1. КОНЦЕПЦИЯ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ

1.1. Технологическая схема имитационного моделирования систем
1.2. Концепция статистического моделирования

1.1. Технологическая схема имитационного моделирования систем

Методологические аспекты и проблемы имитационного моделирования сложных систем удобно рассматривать на укрупненной схеме процесса моделирования (рис.1.1).

Блоки 1,2,3 представляют операции по исследованию одного варианта модели. Эти операции повторяются при различных реализациях случайных процессов, образуя внутренний цикл моделирования (ЦИКЛ I). Процедура выбора оптимального варианта моделируемой системы (блок 4) управляет экспериментом путем изменения соответствующим образом вариантов модели. При этом блоки 1,2,3 (внутренний цикл) охватываются цепями обратной связи (ЦИКЛ II). Связь 3-4-2 отражает адаптацию моделируемой системы. Связь 3-4-1 может возникнуть, если при оптимизации варьируется не только модель системы, но и модель случайных воздействий, рассматриваемых, например, как описание конфликтующей стороны.

Рис. 1.1 Процесс моделирования системы

Оценка результатов исследования вариантов модели оказывается типовой операцией (ЦИКЛ I), многократно выполняемой как в динамическом цикле корректировки модели (ЦИКЛ III), так и в цикле оптимизации (ЦИКЛ II): любой метод поиска экстремума основан на сравнении значений оптимизирующего показателя.

Таким образом, доминирующим в схеме (рис.1.1) является статистическое моделирование (ЦИКЛ 1).

К оглавлению

1.2. Концепция статистического моделирования

В основе статистического моделирования лежит процедура, применяемая для моделирования случайных величин и функций и носящая название метода статистических испытаний (метод Монте-Карло).

Общая схема метода Монте-Карло может быть записана в виде

(1.1)

Результат ищется как математическое ожидание некоторой случайной величины Y, которая чаще всего является неслучайной функцией случайной величины X, имеющей распределение р(х). Нестрогое выражение “случайная величина Х имеет распределение р(х)” и запись Х~ р(х) означают для непрерывной случайной величины, что ее плотность вероятности равна р(х); для дискретной случайной величины функцию р(х) надо понимать как функцию вероятности. Для дискретной случайной величины интеграл (1.1) заменяется суммой S у(х)р(х), в которой суммирование осуществляется по всем возможным значениям Х. Функция у(х) может иметь несколько аргументов, т.е. зависеть от нескольких случайных величин. В таком случае запись (1.1) остается в силе, только интеграл надо считать многомерным, Х рассматривать как вектор, а р(х) – как многомерную плотность (или функцию) вероятности. Приближенная оценка неизвестного математического ожидания, совпадающая с искомым результатом, находится как среднее арифметическое результатов независимых опытов. Это отражено в правой части (1.1). По закону больших чисел среднее арифметическое сходится к математическому ожиданию.

В каждом опыте разыгрывается реализация х случайной величины Х (в i- м опыте реализация xi) в соответствии с распределением р(х) и вычисляется значение функции в виде у(xi). Индекс i подчеркивает, что для каждой (i-й) реализации процесса аргументы, составляющие вектор Х, имеют свои случайные значения. Вычисленное очередное значение у(xi) добавляется к накапливаемой сумме S у(xi). На этом заканчивается очередной опыт. После того как проведено М опытов, вычисляется итоговая оценка в виде правой части выражения (1.1). Опыты повторяются до тех пор, пока дисперсия оценки [image] не снизится до требуемой величины, зависящей от допустимой погрешности и коэффициента доверия.

Проиллюстрируем суть метода Монте-Карло относительно простым примером. Пусть требуется оценить надежность системы (рис.1.2).

Система выполняет свою функцию, если работают цепочки блоков: 1,2,5,7; 1,3,5,7; 1,4,6,7. Какие-то блоки могут отказать. Каждый блок характеризуется временем безотказной работы t i, [image] . Пусть заданы плотности распределения р i(t i), [image] Какова надежность системы в целом?

Рассмотрим случайную величину

g = min { t1, max [min (t4, t6), min [max (t2, t3),]], t7} , (1.2)

где g - время безотказной работы системы. В одном опыте разыгрываются значения всех t i, [image] в соответствии с р i(t i), [image] Используя полученные реализации t i, [image] , по (1.2) вычисляем реализацию g . Один опыт дает одну реализацию (одно выборочное значение) g . Проводим М опытов (испытаний), получаем “статистический” материал (выборку). Берем среднее арифметическое времени безотказной работы системы g с р в качестве оценки надежности системы. При необходимости можно построить закон распределения вероятностей случайной величины g в виде соответствующей гистограммы.

Рис.1.2. Блочная структура системы

Таким образом, испытания реальной системы заменены на испытания математической модели. Каждое испытание сопровождается расчетом. Поэтому имитационное моделирование и называют численным экспериментом на ЭВМ с математической моделью (модель выступает как объект исследования). При реализации испытания возможны и логические операции. И расчетные, и логические операции реализуются на ЭВМ с помощью соответствующих алгоритмов, которые в совокупности и составляют моделирующий алгоритм.

Моделирующий алгоритм обеспечивает построение траекторий смены состояний системы во времени, а воспроизведение случайных факторов, определяющих эти состояния, конструируется с использованием заданных законов случайных событий и величин и реализуется с помощью датчиков базовой случайной величины (БСВ).

К оглавлению

Назад

Содержание

Вперед


2. ПОСТРОЕНИЕ ДАТЧИКОВ БСВ

2.1. Датчики БСВ
2.2. Метод середины квадрата
 2.3. Мультипликативный конгруэнтный метод
2.4. Тестирование равномерности
2.5. Тестирование независимости

2.1. Датчики БСВ

Базовой случайной величиной (БСВ) в статистическом моделировании называют непрерывную случайную величину z, равномерно распределенную на интервале (0,1). Ее плотность распределения вероятностей (п.р.в.) имеет вид:

f (t) = 1, 0 < t < 1. (2.1)

Математическое ожидание (м.о.) и дисперсия БСВ составляют

(2.2)

(2.3)

соответственно.

БСВ моделируется на ЭВМ с помощью датчиков БСВ. Датчик БСВ - это устройство или программа, выдающая по запросу одно или несколько независимых значений z1 , ..., zn БСВ.

Датчики БСВ могут быть трех типов: табличные, физические и программные.

Табличный датчик БСВ - это просто таблица случайных чисел. Основной недостаток такого датчика - ограниченное количество случайных чисел в таблицах. А в статистическом эксперименте часто требуется не ограниченное заранее их количество.

Физический датчик БСВ - это специальное радиоэлектронное устройство в ЭВМ, содержащее источник электронного шума. Шум преобразуется в случайные числа с распределением (1.1). Недостатки физического датчика БСВ: невозможность повторения каких-либо ранее полученных реализаций z1, ... , zn без их предварительной записи в память ЭВМ, схемная нестабильность и сложность тиражирования датчика.

Программный датчик БСВ обычно вычисляет значения z1, z2,..., по какой-либо рекуррентной формуле типа

zi = f ( zn), (2.4)

при заданном стартовом значении z0.

Заданное значение z0 полностью определяет всю последовательность реализаций z1, z2,..., поэтому z часто называют псевдослучайной величиной. Но ее статистические свойства идентичны свойствам "чисто случайной" последовательности, что и обеспечивает успех статистического моделирования.

Программный датчик БСВ имеет следующие преимущества: простота создания датчика, простота применения, простота тиражирования, надежность, быстродействие, высокая точность достижения необходимых статистических свойств, сравнимая с точностью представления вещественных чисел, компактность, повторяемость, когда это нужно, любых последовательностей случайных значений без их предварительного запоминания.

В дальнейшем мы будем рассматривать только программные датчики БСВ.

Имея датчик БСВ z, можно промоделировать любые случайные факторы: непрерывные или дискретные случайные величины (как простые, так и многомерные), случайные события, случайные процессы и поля и т.д. Для этого достаточно соответствующим образом преобразовать последовательность z1, z2, ... . Поэтому БСВ z и называют базовой.

Теоретически в качестве базовой можно было бы взять почти любую случайную величину (с.в.). Использование с.в. z с распределением (2.1) обусловлено технологическими соображениями: простотой и экономичностью датчика, простотой преобразования z в другие случайные факторы, относительной простотой тестирования датчика.

К оглавлению

2.2. Метод середины квадрата

Метод середины квадрата предложен для получения псевдослучайных чисел Д. фон Нейманом в 1946 г. Вот один из вариантов этого метода.

  1. Возьмем произвольное 4-значное число.
  2. Возведем полученное число в квадрат и, если необходимо, добавим к результату слева нули до 8-значного числа.
  3. Возьмем четыре цифры из середины 8-значного в качестве нового случайного 4-значного числа.
  4. Если нужны еще случайные числа, то перейдем к 2.

Например, если взять в качестве начального числа 1994, то из него получается следующая последовательность псевдослучайных чисел: 9760 2576 6357 4114 9249 5440 5936 2360 5696 4444 7491 1150 3225 4006 0480 2304 3084 5110 1121 2566 ...

Сам по себе метод середины квадрата не получил широкого распространения, так как выдает "больше чем нужно малых значений" [26]. Но открытый в нем принцип используется во многих, если не во всех, более поздних датчиках БСВ. Этот принцип состоит в вырезании нескольких цифр из результата какой-либо операции над числами.

К оглавлению

2.3. Мультипликативный конгруэнтный метод

Так называемый мультипликативный конгруэнтный датчик БСВ задается двумя параметрами: модулем m и множителем k. Обычно это достаточно большие целые числа.

При заданных m, k числа z1, z2, ..., вычиcляются по рекуррентной формуле:

Ai = (kAi -1) mod m, i = 1, 2,..., (2.5)

zi = Ai / m,

где m - модуль, k - множитель, A0 - начальное значение, mod - операция вычисления остатка от деления kAi -1 на m.

Таким образом, A1 определяется как остаток от деления k A0 на m; A2 - как остаток от деления kA1 на m и т.д. Поскольку все числа Ai - это остатки от деления на m, то 0 Ј Ai < m. Разделив последнее неравенство на m, видим, что 0 Ј Ai / m< 1, т. е. 0 Ј zi <1.

Из неравенства 0 Ј Ai < m вытекает также, что датчик (2.5) дает периодическую последовательность Ai. Действительно, число всех возможных остатков от 0 до m - 1 равно m и, рано или поздно, на каком-то шаге i обязательно появится значение Ai, уже встречавшееся ранее. С этого момента последовательность Ai “зациклится".

Длина периода T будет не больше m - 1. Например, если встретится остаток Ai= 0, то далее, согласно (2.5), будет Ai+ 1 = 0, Ai+ 2 = 0, ... , т.е. длина периода T = 1. Ненулевых же остатков в интервале 0Ј Ai < m всего m - 1, и, если все они войдут в период, будет T = m - 1. Это имеет место, например, при m = 13, k = 7; в этом случае ряд Ai выглядит так:

1, 7, 10, 5, 9, 11, 12, 6, 3, 8, 4, 2, 1, 7,... .
\_________________________/

T = m - 1 = 12

Поскольку в качестве случайной можно использовать лишь подпоследовательность Ai внутри одного периода, то параметры датчика выбирают так, чтобы длина периода T была максимальной. С учетом ограничения TЈ m - 1 модуль m берут максимально возможным. Чтобы упростить вычисление остатков по (2.5), для двоичных ЭВМ часто берут m = 2n. Рекомендуется также брать достаточно большой множитель k, причем взаимно простой с m.

В [30] можно найти подробные рекомендации по выбору параметров m, k и начального значения A0 . Заметим, однако, что в настоящее время не известны правила, которые гарантировали бы высокое качество датчика без его специального статистического тестирования.

Датчик (2.5) называют мультипликативно-конгруэнтным потому, что он использует две основные операции - умножение (англ. multiplication) и вычисление остатка (в теории чисел - получение конгруэнтного числа). Можно было бы поэтому перевести его название и как "множительно-остатковый датчик".

Обратим внимание также и на то, что операция вычисления остатка воплощает здесь упоминавшийся в п. 2.2 неймановский принцип вытаскивания цифр. Это становится очевидным, если записывать числа в системе счисления с основанием m. Тогда операция X mod m означает выбор последней цифры из числа X. Для m = 2n операция X mod m означает также выделение последних n цифр из двоичной записи числа X.

В заключение параграфа приведем таблицу параметров датчиков, предлагаемых в некоторых публикациях и программных продуктах.

Таблица 2.1

Место использования датчика (программный продукт или публикация) | модуль m | множитель k
Язык моделирования СИМУЛА | 235 | 516
Библиотека АЛГОЛ-процедур | 235 | 515
Пакеты LLRANDOM, IMSL | 231 - 1 (простое число) | 16807
Язык моделирования SIMSCRIPT | 231 - 1 | 63036001
[26, 5] | 240 | 517
[26] | 236 | 515
[27] | 242 | 517

К оглавлению

2.4. Тестирование равномерности

Обозначим равномерное распределение вероятностей на интервале (0,1) через R[0,1]. Тогда утверждение, что БСВ z имеет распределение R[0,1], можно кратко записать в виде z ~ R[0,1].

С помощью статистических тестов проверяют два свойства датчика, делающих его точной моделью идеальной БСВ, - это равномерность распределения чисел zi, выдаваемых датчиком на интервале (0,1), и их статистическая независимость. При этом числа zi рассматривают как реализации некоторой с.в., т.е. как статистическую выборку.

Достаточно простым методом проверки равномерности распределения является частотный тест. Он основан на законе больших чисел и выполняется по следующему алгоритму.

1. Разобьем интервал (0,1) на K равных отрезков (например, K = 10).

2. Сгенерируем n чисел z1,..., zn с помощью тестируемого датчика БСВ (например, n = 100).

3. Подсчитаем, сколько чисел попало в каждый из k отрезков, т.е. найдем числа попаданий n1,...,nk.

4. Рассчитаем относительные частоты попаданий в отрезки: ,..., [image].

5. Построим гистограмму частот [image] 1,..., [image] k на K отрезках интервала (0,1).

6. Повторим действия (2) - (5) для большего значения n (например, для n =10 000).

7. Оценим по полученным гистограммам сходимость каждой частоты [image] i к вероятности p = 1/K того, что БСВ попадет в i-й отрезок. Согласно закону больших чисел должно быть

n ® Ґ (2.6)

Это значит, что высоты столбиков во второй гистограмме должны в целом быть ближе к уровню 1/K, чем в первой.

Тестирование датчика на равномерность можно совместить с оцениванием м.о. и дисперсии с.в. Оценки [image] и [image] для м.о. и дисперсии рассчитываются соответственно по формулам:

(2.7)

(2.8)

С ростом n оценки [image] и [image] должны сходиться по вероятности к точным значениям M(z) = 1/2, D(z) = 1/12 = 0.08333... .

В табл. 2.2 приводится пример тестирования одного из распространенных датчиков БСВ.

Таблица 2.2

Пример испытания датчика
Рассчитываемые оценки |  | Длина выборки n
 |  | n =1000 | n = 100 000
Оценки м.о. и дисперсии |  | 0.4928 | 0.4992
 |  | 0.0841 | 0.0832
Частоты | 1 | 0.114 | 0.10136
 | 2 | 0.097 | 0.09974
 | 3 | 0.100 | 0.09868
 | 4 | 0.108 | 0.10130
 | 5 | 0.094 | 0.10155
 | 6 | 0.096 | 0.09896
 | 7 | 0.096 | 0.09928
 | 8 | 0.102 | 0.10015
 | 9 | 0.104 | 0.10085
 | 10 | 0.089 | 0.09813
Сумма частот |  | 1.000 | 1.000

На рис. 2.1 и 2.2 показаны соответствующие гистограммы. Судя по этим результатам, протестированный здесь датчик отвечает требованию равномерности.

Рис. 2.1. Гистограмма для n=1000.

Рис.2.2. Гистограмма для n =100000

К оглавлению

2.5. Тестирование независимости

Простейшую проверку статистической независимости реализаций z1, z2, ..., можно осуществить, оценивая корреляцию между числами zi и zi+ s, отстоящими друг от друга на шаг s >1.

Для вывода формулы, по которой можно рассчитать коэффициент корреляции чисел zi и zi+ s , рассмотрим две произвольные с.в. x, y. Коэффициент корреляции определяется для них формулой:

(2.9)

Если известно, что x, y ~ R[0,1], то M(x) = M(y) = 1/2 и D(x) = D(y) =1/12, то есть (2.9) принимает вид:

R(x,y) = 12 M(xy) - 3. (2.10)

Условимся рассматривать пару чисел (zi , zi+ s) как реализацию пары с.в. (x,y). Тогда в выборке z1,..., zn имеем всего n - s реализаций этой пары:

(z1 , z1+ s),(z2 , z2+ s),... ,(zn , zn+ s).

По ним можно рассчитать оценку R' коэффициента корреляции R(x,y), заменяя в (1.10) м.о. M(xy) соответствующим средним арифметическим:

(2.11)

С ростом n оценка R' должна приближаться к нулю, в противном случае датчик БСВ не отвечает требованию независимости.

Конечно, если R' сходится к нулю, то это еще не гарантирует наличие независимости, но все же один из тестов оказывается успешно выдержанным. При желании всегда можно продолжить испытания датчика другими методами.

В табл. 2.3 показан пример тестирования датчика БСВ. Как видно из таблицы, оценки коэффициентов корреляции для s =1,2,3 постепенно сходятся к нулю. На рис. 2.3 изображены соответствующие графики.

Таблица 2.3

n | Оценка R' для s = 1,2,3
 | s = 1 | s = 2 | s = 3
1000 | 0.175 | 0.142 | 0.232
2000 | 0.071 | 0.059 | 0.111
3000 | 0.007 | 0.021 | 0.043
4000 | 0.017 | 0.040 | 0.043
5000 | 0.019 | 0.045 | 0.042
6000 | 0.009 | 0.031 | 0.039
7000 | -0.003 | 0.019 | 0.027
8000 | -0.003 | 0.021 | 0.028

Рис.2.3. Уточнение коэффициентов корреляции

В настоящее время существуют десятки специальных методов, пригодных для тестирования датчиков БСВ. Однако эта проблема имеет не математическую природу, поэтому она не имеет корректного математического решения [9].

Дж. Клейнен в [12] отмечает, что "Наилучший в настоящее время, по-видимому, мультипликативный генератор, однако он требует тщательного выбора сомножителя".

Т.Дж. Шрайбер в [29] описал датчик БСВ, встроенный в язык GPSS (в версиях для IBM 360/370).

Точное описание алгоритма и данные тестирования этого датчика в [29] не приводятся, как и в фирменной документации [6, 24]. Наше доверие к датчику может питаться авторитетом фирмы IBM и, в особенности, широким апробированием языка GPSS в прикладных исследованиях и в учебных курсах моделирования.

Хорошим подтверждением высокого качества датчика можно считать также работу [5], где на GPSS моделируются системы массового обслуживания и результаты сравниваются с точными аналитическими решениями.

К оглавлению


3. МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ФАКТОРОВ

3.1. Случайные факторы и их модели
3.2. Моделирование случайных событий
3.3. Моделирование дискретных случайных величин
3.4. Моделирование непрерывных с.в. методом обращения
3.5. Моделирование экспоненциальной с.в.
3.6. Доказательство метода обращения
3.7. Линейные преобразования с.в.
3.8. Моделирование равномерной с.в.
3.9. Моделирование эрланговской с.в.
3.10. Моделирование нормальной с.в.

3.1. Случайные факторы и их модели

Под случайным фактором условимся понимать случайное событие, случайную величину, многомерную случайную величину, случайную функцию одного или нескольких переменных и вообще любые объекты, выбор которых определяется соответствующими вероятностными распределениями. А под реализацией случайного фактора - сам акт выбора, т.е. статистическую выборку одного элемента из множества, наделенного вероятностной мерой.

В данной главе рассмотрение будет ограничено моделированием случайных событий и случайных величин. Технология моделирования более сложных случайных факторов может быть изучена по [30, 5].

При моделировании случайных факторов будем исходить из того, что используемый датчик БСВ является идеальной моделью базовой случайной величины z ~ R[0,1], т.е. выдает равномерно распределенные и независимые случайные значения z i. Это позволит нам строить модели случайных факторов и устанавливать их свойства теоретически - опираясь только на положения теории вероятностей и не прибегая к громоздким процедурам статистического тестирования.

В то же время в целях наглядности модели случайных факторов будем иллюстрировать примерами с простой статистической обработкой выходных данных.

К оглавлению

3.2. Моделирование случайных событий

Пусть задана некоторая система непересекающихся событий {A 1,...,A n} c их вероятностями p 1,..., p n соответственно; p 1+...+ p n= 1.Чтобы построить программную модель такого события, разобьем интервал (0,1) на n отрезков длиной p 1, p 2, ..., p n. Это всегда возможно, так как p 1+...+ p n= 1. Тогда модель случайного события A может быть написана в соответствии с простым алгоритмом, изображенным на рис. 3.1, где запись (p i) обозначает "отрезок длиной p i".

В качестве примера построим модель операции, состоящей в вытаскивании шара из урны, содержащей пять белых шаров (Б), три красных (К) и два черных (Ч). Так как исходы Б,К,Ч имеют вероятности p 1 = 0.5, p 2 = 0.3 и p 3 = 0.2 соответственно, то интервал (0,1) разбиваем на отрезки (0;0.5), (0.5;0.8) и (0.8;1).

Алгоритм моделирования имеет следующий вид:

(1) Получить значение z из датчика БСВ.
(2) Если z Ј 1/2 , вывести "Б", иначе, если z Ј 8/10 , вывести "К", иначе вывести "Ч".

Рис. 3.1. Модель случайного события

Вот пример 60-кратного выполнения этого алгоритма на компьютере; мы видим, что частота появления каждого исхода примерно соответствует его вероятности:

БКБКБКБКББКЧБББККББКЧЧББЧЧКББЧБЧББЧБКЧЧБББККББЧКБЧКББЧКБЧББ.

Так, исход "Б" здесь появился 31 раз (52% случаев), "К" - 15 раз (25%) и "Ч" - 14 раз (23%).

К оглавлению

3.3. Моделирование дискретных случайных величин

Моделирование дискретной с.в. в принципе ничем не отличается от моделирования случайного события. Дискретная с.в. d задается множеством возможных значений {d 1, ... , d n} и их вероятностями p 1, ... , p. Поэтому дискретная с.в. реализуется по тому же алгоритму рис.2.1, только вместо "Исход A i" в нем нужно записать "d = d i" для всех i = 1, ... , n.

Единственное, что здесь имеет смысл добавить сверх сказанного, это о возможности модифицировать алгоритм рис. 3.1 в двух частных случаях в целях упрощения программы.

Первый случай - это когда dО {0, 1, ... , n - 1} и исходы равновероятны: p 1= ... = p n = 1/n. В этом случае все промежуточные операции, включая разбиение интервала (0,1) на n отрезков и поиск среди них того, в который попало значение z, сокращаются за счет применения функции л x ы вычисления целой части от x. И весь алгоритм сворачивается в одну формулу:

d = л nz ы . (3.1)

Действительно, в этом случае при zО (0,1/n) будет d = 0, при zО (1/n,2/n) будет d = 1 и т.д. Случаи, когда z попадает точно на границы отрезков, здесь не рассматриваются, так как их вероятность равна нулю.

Ясно, что (3.1) легко приспособить к близким ситуациям. Например, для моделирования игральной кости с числом очков от 1 до 6 можно ее "выбрасывание" реализовать по формуле

d = л 6z ы + 1, (3.2)

где z получается от датчика БСВ, d - выпавшее число очков.

Применение (3.2) для 60 "выбрасываний кости" компьютером дало последовательность:

321236246263321315464215224214664122142366236246611664533544.

Здесь цифры 1,2,3,4,5,6 выпали соответственно 9,14,9,11,4 и 13 раз, т.е. в процентах это составило 15%, 23%, 15%, 18%, 7% и 22%.

Второй случай, когда алгоритм рис. 3.1 следует упрощать, касается дискретной с.в. с бесконечным множеством возможных исходов. В такой ситуации часто можно получать компактную программу, применяя рекурсивный (циклический) вариант алгоритма.

В качестве примера построим генератор пуассоновской с.в. Пуассоновская с.в. может принимать значения {0, 1, ... , K, ...}, вероятности которых определяются по формуле

, (3.3)

где a - параметр распределения, a > 0.

Убедимся, что сумма вероятностей исходов равна единице:

. (3.4)

Разобьем мысленно интервал (0,1) на бесконечное число отрезков длиной p0 , p1 , ... , p k, ... . Согласно рис. 3.1 мы должны обратиться к датчику БСВ и для полученного значения z найти отрезок, в который оно попало. Проверить для первого отрезка условие zО ( p 0) можно с помощью неравенства zЈ Г0, где Г0 - правая граница интервала (p 0), т.е. Г0= р0= е-а. Если после отрезка (p k) нужно проверять отрезок (p k+1), то для этого используем неравенство z Ј Гk+1, где Гk+1 -правая граница отрезка (p k+1) - получается по формуле Гk+1= Гk + p k+1. Вычисление p k+1 будем выполнять не непосредственно по (3.3), а по соответствующему рекуррентному выражению

, (3.5)

чтобы не считать факториалов и степеней.

Алгоритм, соответствующий высказанным соображениям, представлен на рис. 3.2. Переменные Г и р, изменяющиеся при повторении цикла, соответствуют величинам Гk+1 и p k+1. Полученный алгоритм имеет достаточно компактный вид. Что касается быстродействия алгоритма, то его можно охарактеризовать числом N повторений цикла. Очевидно, что N = d + 1, где d - вычисленное случайное значение пуассоновской с.в. Поскольку M(d) = a, то среднее число повторений цикла составляет

M(N) = M(d + 1) = a + 1. (3.6)

Нетрудно определить и дисперсию D(N):

D(N) = D(d + 1) = D(d) = a. (3.7)

Рис.3.2. Схема генератора пуассоновской с.в.

На рис. 3.3 приведена частотная гистограмма для 100 000 реализаций пуассоновской с.в., полученных с помощью описанного генератора.

Рис. 3.3. Пуассоновская с.в. при а=3..5

К оглавлению

3.4. Моделирование непрерывных с.в. методом обращения

Пусть требуется реализовать непрерывную с.в. x, которая имела бы заданную функцию распределения вероятностей (ф.р.в.) F(t) = P{x Ј t}, т.е. требуется, чтобы было x ~ F(t). Согласно методу обращения это можно сделать с помощью БСВ z по следующей формуле:

x = F-1(z) . (3.8)

Примечание. Выражение x = F-1(z) для расчета x через z можно получить следующим известным способом:

(1) Записать формальное уравнение F(x) =z.
(2) Разрешить его относительно x.

К оглавлению

3.5. Моделирование экспоненциальной с.в.

Экспоненциальная с.в. x имеет ф.р.в.

F(t) = 1 – e-l t , (3.9)

где t і 0 и параметр l > 0 . Для экспоненциальной с.в. x

M(x) = 1 / l , D(x) = 1 / l 2. (3.10)

Для построения генератора такой с.в. используем метод обращения.

  1. Записываем формальное уравнение F(x) = z: 1 – e-l t = z.
  2. Решаем его относительно x:

x = - (1/l ) ln (1- z). (3.11)

Формулу (3.11) можно упростить, заменив (1 - z) на z, так как обе эти величины совпадают по распределению. Тогда из (3.11) получаем

x = - (1 / l ) ln (z). (3.12)

Рис. 3. 4. Частотная гистограмма
экспоненциальной с.в. (м.о. = 1) при n = 10000

При независимых z генератор (3.12) дает независимые значения экспоненциальной с.в. x. На рис. 2.4 приводится частотная гистограмма, полученная при испытании этого генератора для l = 1. Оценки моментов составили [image] = 0.990, [image] = 0.978.

К оглавлению

3.6. Доказательство метода обращения

Требуется доказать, что если x = F -1(z), где F - ф.р.в., то x ~ F(t). Предположим, что в действительности x ~ F1(t). Равенство F-1 и F доказывается следующей цепочкой эквивалентных выражений:

F1 (t) = Р{x Ј t} = P{F-1(z) Ј t} = P{F(F-1 (z))Ј F(t)} = P{z Ј F(t)} = P{z О(0,F(t))} = F(t).

Первое равенство здесь имеет место по определению, последнее - в силу равномерного распределения z на отрезке (0,1), на котором лежит и отрезок (0,F(t)) длиной F(t). Остальные равенства - по причине эквивалентности соответствующих событий, помещенных под знак вероятности (эти события суть одно и то же неравенство, переписываемое в разной форме).

К оглавлению

3.7. Линейные преобразования с.в.

Пусть x - с.в. и

y = Ax + B, (3.13)

где A, B - константы; A - коэффициент, B - смещение. Преобразование с.в. x вида (3.13) называется линейным.

Многие непрерывные с.в. в результате линейного преобразования не меняют вид распределения (меняются только его параметры). К таким относятся экспоненциальная или эрланговская с.в., вид распределения которой сохраняется только при масштабном преобразовании, т.е. когда смещение B = 0.

Очевидно, при линейном преобразовании (3.13) имеет место:

M(y) = AM(x) + B, (3.14)

D(y) = A2D(x). (3.15)

Линейное преобразование (3.13) часто применяется для того, чтобы при наличии генератора непрерывной с.в. x получить другую с.в. y с тем же распределением, но с другими значениями м.о. и дисперсии.

К оглавлению

3.8. Моделирование равномерной с.в.

Равномерно распределенную с.в. x можно реализовать путем линейного преобразования БСВ z:

x = Az + B. (3.16)

При этом, поскольку z ~ R[0,1], а (2.16) увеличивает диапазон значений с.в. в A раз и смещает вправо на величину B, то x будет иметь диапазон значений (B, B + A), т.е.

x ~ R[B, B + A]. (3.17)

Кроме того, из M(z) = 1/2, D(z) = 1/12 и из (2.14),(2.15) вытекает:

M(x) = A/2, (3.18)

D(x) = A2/12. (3.19)

К оглавлению

3.9. Моделирование эрланговской с.в.

Эрланговская с.в. имеет ф.р.в.

, t > 0, (3.20)

где k - порядок распределения.

В данном случае не удается найти простое выражение для обратной функции F-1, поэтому применить метод обращения так, как он описан в п. 2.4, невозможно. В подобных ситуациях для построения генератора с.в. могут пригодиться следующие методы:

а) приближенное обращение функции F(t) путем ее разложения в ряд или путем аппроксимации другими функциями;

б) приближенное численное обращение функции F(t) по таблице ее значений с применением линейной интерполяции (этот метод реализован в GPSS);

в) точная реализация с.в. x, как функции от других с.в. (не обязательно от базовой с.в., как в методе обращения);

г) приближенная реализация x как функции от других с.в.;

д) точная или приближенная реализация x по функции f(t) = dF(t)/dt - т.е. по плотности распределения вероятностей (п.р.в.).

В данном случае для реализации эрланговской с.в. используем вариант в).

Итак, известно, что случайной величиной Эрланга k-го порядка называется сумма k независимых экспоненциальных с.в., имеющих одно и то же значение параметра l . Согласно этому определению и с учетом (3.12) генератор эрланговской с.в. x будет просто вычислять ее как сумму:

(3.21)

где zi - реализации БСВ, i = 1, ...,K.

На рис. 3.5 приводятся результаты вычисления 10 000 значений x по (2.21) при l = 1, K = 3. Так как здесь x = x1 + x2 + x3 , где M(xi) = 1/l =1, D(xi)=1/l 2 = 1, i = 1,2,3, то M(x)=3, D(x)=3. Статистические оценки, полученные при испытаниях, составили [image] (x) = 2.968, [image] (x) = 2.958.

Рис. 3.5. Гистограмма эрланговской с.в. 3-го порядка при n = 10000

К оглавлению

3.10. Моделирование нормальной с.в.

Утверждение, что некоторая с.в. x имеет нормальное распределение с м.о. M(x) = m и дисперсией D(x) = s 2, будем записывать в виде x ~ N [m, s ]. Плотность распределения вероятностей (п.р.в.) такой с.в. имеет вид

(3.22)

Применение метода обращения для моделирования нормальной с.в. затрудняется в связи с тем, что ее ф.р.в. F(t) не выражается в элементарных функциях в виде конечной формулы. Поэтому здесь часто используют упомянутый в п. 3.9 вариант г). А именно, нормальную с.в. часто реализуют приближенно как сумму других с.в., основываясь на центральной предельной теореме теории вероятностей.

Рассмотрим сумму S = x1 +...+ xn произвольных с.в., которые имеют одинаковое распределение с конечными м.о. M и дисперсией D. Яcно, что M(S) = nM и D(S) = nD. Поэтому для линейного преобразования (нормированного значения) суммы (S - nM) [image] имеем

, (3.23)

. (3.24)

Согласно центральной предельной теореме при x ® Ґ с.в. (S-nM)/ [image] стремится по распределению к нормальной с.в. x ~ N[0,1]. Поэтому приближенный генератор с.в. x ~ N[0,1] можно построить по формуле x(S-nM)/ [image] , взяв достаточно большое n. Хорошая точность получается, например, если взять сумму S двенадцати слагаемых z1 + ... + z12, в качестве которых использовать независимые реализации БСВ. Тогда M = M(z) = 1/2, D = D(z) = 1/12 и x = (S- 12M)/ [image] = S- 6, т.е.

(3.25)

На рис. 3.6. показана частотная диаграмма, которая получена для 100 000 значений x, вычисленных по формуле (3.25). При этом статистические оценки, полученные для м.о. и дисперсии, составили соответственно [image] = - 0.0085 и [image] =1.0029.

Используя линейное преобразование с.в. x, можно получить нормальную с.в. y с любыми значениями m, s .

Рис. 3. 6. Гистограмма нормальной с.в. при

n = 100000

К оглавлению


4. ПЛАНИРОВАНИЕ СТАТИСТИЧЕСКОГО ЭКСПЕРИМЕНТА

4.1. Задача планирования
4.2. Теоретическое решение задачи планирования
4.3. Метод автоостанова
4.4. Интерактивный контроль

4.1. Задача планирования

Пусть X = (x1, ... , xk ) - вектор случайных величин, имеющий распределение PX . Во многих случаях цель статистического эксперимента может быть сведена к определению математического ожидания с.в. y = g(X), где g(X) - заданная функция с.в. X. Поэтому типовая схема статистического эксперимента имеет вид рис. 4.1.

Рис. 4.1. Типовая схема статистического эксперимента

Здесь оценка [image] приближается к точному значению M(y) с ростом числа опытов n.

Задача планирования статистического эксперимента состоит в нахождении такого n, при котором достигается заданная точность оценки [image] .

К оглавлению

4.2. Теоретическое решение задачи планирования

Оценка [image] является случайной величиной, так как представляет собой функцию от случайных величин:

. (4.1)

В типовой схеме статистического эксперимента (рис. 4.1) используются независимые реализации xi , поэтому с.в. yi в (4.1) также независимы. Из центральной предельной теоремы вытекает, что с.в. в (4.1) при больших n имеет нормальное распределение, т.е. ~N [m,s ]. Определим параметры этого распределения:

, (4.2)

, (4.3)

. (4.4)

Таким образом, при увеличении числа опытов на порядок s уменьшается в [image] » 3 раза. На рис. 4.2. схематически показан вид распределения оценки [image] в зависимости от n.

Рис. 4.2. Распределение с.в. в зависимости от n

Диапазон вероятных отклонений оценки от точного значения M(y) сужается пропорционально [image]. Параметр s используют как показатель точности оценки. Поскольку [image] имеет нормальное распределение, то практически достоверно, что [image] отклоняется от искомого M(y) не более, чем на 3s . Можно сказать, что s является аналогом абсолютной погрешности, а 3s - самой абсолютной погрешностью.

В качестве аналога относительной погрешности для с.в. можно рассматривать ее коэффициент вариации:

v = [image], (4.5)

а в качестве самой относительной погрешности - величину 3v. Из (4.2) и (4.4) находим, что

v = [image] = [image], (4.6)

где vy - коэффициент вариации с.в. y.

Соотношения (4.4) и (4.6) показывают, что "абсолютная погрешность" 3s и "относительная погрешность" 3v оценки [image] убывают пропорционально [image] .

Из (4.4) находим решение задачи планирования: для достижения заданной точности s требуется

(4.7)

опытов. Если требование к точности задано в форме коэффициента вариации v, то требуемое число опытов определяется из (4.6) в виде:

. (4.8)

В решении (4.7) кроме заданного s для определения n необходимо знать s y, в варианте (4.8) - коэффициент вариации vy. Ни то, ни другое до эксперимента, как правило, не известно. Поэтому планирование числа опытов на практике осуществляется в ходе самого статистического эксперимента. Достаточно удобными методами такого планирования являются метод автоостанова и метод интерактивного контроля.

К оглавлению

4.3. Метод автоостанова

Метод автоостанова покажем на примере достижения точности v = 0.01. Схема статистического эксперимента с автоостановом изображена на рис. 4.3.

Как видно из рисунка, в ходе эксперимента ведется непрерывный контроль за оценкой v' коэффициента вариации v. Когда оценка v' устойчиво входит в зону v' < 0.01, эксперимент завершается.

Устойчивое достижение неравенства v'< 0.01 здесь обеспечивается требованием, чтобы это неравенство было подтверждено A = 10 раз.

Заметим, что минимальное начальное значение счетчика подтверждений на рис. 4.3, задаваемое в начале алгоритма, должно составлять A = 2, так как необходимо исключить одно ложное подтверждение, которое имеет место при n = 1. Действительно, при n = 1 оценка дисперсии равна нулю, поэтому получается v' = 0. При A = 10 получится несколько "лишних" повторений цикла. Однако, так как n обычно составляет десятки тысяч, увеличение его на несколько лишних единиц не имеет существенного практического значения.

Рис. 4.3. Схема эксперимента с автоостановом

К оглавлению

4.4. Интерактивный контроль

В случае, когда длительность статистических экспериментов не слишком велика, удобен также вариант планирования, который состоит в простом визуальном наблюдении за изменением точности оценок в ходе эксперимента. Значения s ' или v' выводятся для некоторых n на экран компьютера и по ним принимается решение об остановке или продолжении эксперимента


5. СТАТИСТИЧЕСКИЙ РАСЧЕТ НАДЕЖНОСТИ

5.1. Классическая задача надежности
5.2. Аналитический метод расчета
5.3. Статистический метод расчета
5.4. Пример статистического расчета
5.5. Сравнение аналитического и статистического методов

5.1. Классическая задача надежности

В классической задаче надежности рассматриваются системы следующего вида.

Система состоит из n элементов. Каждый элемент случайно может находиться либо в рабочем состоянии, либо в состоянии отказа. Совокупность состояний элементов удобно выражать двоичным вектором X:

X = (x1, x2, ... , xn ), (5.1)

где xi - состояние i-го элемента, xi О{0,1}. При этом xi = 1 означает, что i-й элемент находится в состоянии отказа, xi= 0 - что он работает. С.в. xi независимы, так как элементы отказывают независимо друг от друга. Известны вероятности отказа элементов: pi = P{xi=1} = M(xi). Обозначим набор вероятностей через P:

P = ( p1, p2 , ... , pn ) (5.2)

Так как вектор P полностью задает закон распределения с.в. X, то будем пиcать X ~ P.

Система в целом также либо работает, либо находится в состоянии отказа. Состояние системы однозначно определяется совокупностью состояний элементов, т.е. является логической функцией вида

y = y(x1, x2 , ... , xn ) = y(X), (5.3)

где y О{0,1} - состояние системы, y = 1 - означает состояние отказа, y = 0 - рабочее состояние системы.

Обозначим через Q вероятность отказа системы, т.е. Q = P{y=1}= M(y).

В классической задаче надежности требуется найти вероятность отказа системы Q при известных вероятностях отказа элементов (5.2).Таким образом, чисто формально условие классической задачи надежности может быть изложено следующим образом.

Даны n-мерная двоичная с.в. X ~ P и логическая функция y = y(X). Требуется найти значение Q = M(y).

Как видим, данная задача вполне может рассматриваться как типовая задача статистического эксперимента.

К оглавлению

5.2. Аналитический метод расчета

5.2.1. Расчет приводимых схем

Аналитическое решение классической задачи надежности в значительной степени опирается на представление зависимости y(X) в виде надежностного графа. Например, на рис. 5.1 представлен граф некоторой системы, состоящей из 4 элементов. Считается, что отказу i-го элемента соответствует пропадание i-го ребра графа (ребра изображаются в виде прямоугольников). Система в целом считается отказавшей в тех случаях, когда между полюсами A, B не остается ни одного пути.

Такой граф определяет состояние системы для любых состояний элементов, т.е. полностью определяет функцию y(X).

Если все элементы в графе параллельны (рис. 5.2), то система отказывает тогда и только тогда, когда отказывают все элементы (все xi = 1), т.е.

y = [image]. (5.4)

Следовательно, для параллельного соединения

Q = [image]. (5.5)

Здесь Х - символ произведения, а [image] в силу независимости xi, i = 1...n.

Рис. 5.1. Пример надежностного графа Рис. 5.2. Параллельное соединение элементов

Если элементы соединены последовательно (рис.5.3), то система работает (y = 0) тогда и только тогда, когда работают все элементы (все xi = 0). Этот факт можно выразить формулой

y = 1 - [image] (5.6)

Следовательно, для такого соединения

Q = M(y) = 1 - ХM(1 - x i) = 1 - Х(1 - p i) = 1 - Хq i, (5.7)

где q i = 1 - pi - вероятность работоспособного состояния i-го элемента, i = 1...n.

Рис. 5.3. Последовательное соединение

Используя (5.5) и (5.7), можно рассчитать надежность любой системы, изображенной в виде параллельно-последовательного (т.е. так называемого приводимого) графа.

Например, для рис 5.1 "параллельная" формула (5.5) применяется к параллельному соединению ребра 4 и группы ребер 1, 2, 3. При этом вероятность отказа для групп 1, 2, 3 выражается с помощью "последовательной" формулы через p1 и через вероятность отказа группы 2, 3. А последняя вероятность - прямо через вероятности p2 p3:

Q =M(y) = p4p123
= p4(1– q1q23) =
p4(1– q1(1- p23)) = p4(1–
q1(1- p2p3)).             (5.8)

Эта формула, полученная для рис. 5.1, содержит четыре вхождения аргументов. Будем говорить, что ее длина равна 5. Длина расчетной формулы для приводимого графа из n элементов составит n. Как видим, формулы приводимых графов получаются весьма компактными и простыми для вычислений.

Кроме того, процесс получения расчетной формулы здесь также весьма прост. Он не требует анализа функции y(X) или путей A,B. Вместо этого выполняется разбиение графа на подграфы.

5.2.2. Неприводимые схемы

Однако аналитический метод расчета значительно усложняется, если надежностный граф оказывается неприводимым (как на рис. 5.4), или если он оказывается многополюсным (как на рис. 5.5), или если зависимость y(X) вообще не может быть изображена надежностным графом.

В подобных случаях для получения точной формулы расчета M(y) часто приходится анализировать все множество элементарных исходов с.в. X= (x1,..., xn) с тем, чтобы определить, какие из них дают результат y(X) = 1, и найти для таковых сумму их вероятностей, составляющую значение Q = M(y). Например, для рис. 5.4 имеем 25 = 32 исхода (см. табл. 5.1). Вероятность отказа системы определяется как сумма вероятностей тех исходов, в которых y = 1.

Так, исход X = 01001 означает, что x1= 0, x2= 1, x3= 0, x4= 0, x5 = 1 (т.е. элементы 2 и 5 отказали). Смотрим на рис. 5.5. Мы видим, что при удаленных ребрах 2 и 5 полюса несвязны. Значит, для данного исхода y = 1.

Суммируя вероятности всех исходов с y = 1, получаем для рис. 5.4 следующую расчетную формулу:

M(y) = q1p2q3q4p5 + q1p2q3p4p5 + q1p2p3q4p5 + q1p2p3p4p5 + p1q22q3p4q5 + p1q2q3p4p5 + p1q2p3q4p5 + p1q2p3p4q5 + p1q2p3p4p5 + p1p2q3q4p5 + p1p2q3p4q5 + p1p2q3p4p5 + p1p2p3q4p5 + p1p2p3p4q5 + p1p2p3p4p5

или после упрощения:

Q = M(y) = p1p4+ p2p5(1- p1p4) + p3 (q1 p2 p4q5+ p1q2q4p5). (5.9)

Рис.5.4. Пример неприводимого графа

Таблица 5. 1

Значения с.в. X= (x1x2x3x4x5) | y(X) | Вероятности исходов для y = 1
0 0 0 0 0 | 0 |
0 0 0 0 1 | 0 |
0 0 0 1 0 | 0 |
0 0 0 1 1 | 0 |
0 0 1 0 0 | 0 |
0 0 1 0 1 | 0 |
0 0 1 1 0 | 0 |
0 0 1 1 1 | 0 |
0 1 0 0 0 | 0 |
0 1 0 0 1 | 1 | q1p2q3q4 p5
0 1 0 1 0 | 0 |
0 1 0 1 1 | 1 | q1p2q3p4p5
0 1 1 0 0 | 0 |
0 1 1 0 1 | 1 | q1p2p3q4p5
0 1 1 1 0 | 1 | q1p2p3p4q5
0 1 1 1 1 | 1 | q1p2p3p4p5
1 0 0 0 0 | 0 |
1 0 0 0 1 | 0 |
1 0 0 1 0 | 1 | p1q2q3p4q5
1 0 0 1 1 | 1 | p1q2q3p4p5
1 0 1 0 0 | 0 |
1 0 1 1 0 | 1 | p1q2p3p4q5
1 0 1 1 1 | 1 | p1q2p3p4p5
1 1 0 0 0 | 0 |
1 1 0 0 1 | 1 | p1p2q3q4p5
1 1 0 1 0 | 1 | p1p2q3p4q5
1 1 0 1 1 | 1 | p1p2q3p4p5
1 1 1 0 0 | 0 |
1 1 1 0 1 | 1 | p1p2p3q4p5
1 1 1 1 0 | 1 | p1p2p3p4q5
1 1 1 1 1 | 1 | p1p2p3p4p5

Самой сложной частью аналитического расчета надежности в подобных случаях является определение значения y для всех 2 n значений X, ведь для этого нужно, вообще говоря, 2 n раз проанализировать связность полюсов графа. Или, если зависимость y(X) задана не графом, а, например, некоторым алгоритмом, то нужно 2 n раз выполнить этот алгоритм.

5.2.3. Многополюсный граф

На рис. 5.5 приведен многополюсный надежностный граф, изображающий бортовую систему электрообеспечения [25].

Аналитический расчет этой довольно простой на вид системы в литературе отсутствует (см., например, [16,25]) в силу того, что он требует анализа 215= 32 768 исходов.

Очевидно, если неприводимая система содержит не 15, а несколько десятков или сотен элементов, то ее аналитический расчет может стать практически невозможным

.

P1 = P2 = P3 = 10-5
P5 = P7 = P8 = 3 * 10-5
P4 = P6 = P9 = 2 * 10-5
P10 = ... = P15 = 5 * 10-5

      Рис. 5.5. Граф системы электрообеспечения

К оглавлению

5.3. Статистический метод расчета

Статистический расчет M(y) выполняется по типовой схеме статистического эксперимента (рис. 5.6). При этом с.в. X является n-мерной двоичной, а функция y(X) вычисляется по графу или по любой другой форме задания зависимости y(X). Время вычисления y(X) в одном опыте по графу примерно пропорционально числу элементов системы.

Рис. 5.6. Схема статистического расчета надежности

В каждом из N опытов значение с.в. X = (x1,...,x2) разыгрывается по следующей схеме:

    IF RANDOM < p1 THEN x1:=1 ELSE x1:=0;

    IF RANDOM < p2 THEN x2:=1 ELSE x2:=0;

    .

    .

    .

    IF RANDOM < pn THEN xn:=1 ELSE xn:=0;

где RANDOM – независимое значение БСВ. Или, если вероятности p1, p2, ... , pn предварительно записаны в массив P, то x1, x2, ... , xn разыгрываются с помощью цикла:

    FOR i = 1 TO n DO

    IF RANDOM < P[i] THEN x[i]:=1 ELSE x[i]:=0.

Оценим затраты машинного времени на весь расчет надежности. Пусть, например, требуется получить оценку [image] вероятности Q = M(y) с точностью v = 0.1, где v – коэффициент вариации оценки. Обычно такая точность достаточна для практических целей. Тогда, согласно (3.8), число опытов

(5.10)

где v2y = D(y)/M2 (y)=Q(1-Q)/Q2=(1-Q)/Q. Как правило, вероятность отказа системы Q мала, поэтому v2y= (1- Q)/Q » 1/Q, и (5.10) принимает вид:

. (5.11)

Здесь мы видим замечательное свойство статистического метода расчета: число опытов N зависит от надежности Q и точности v, но не зависит от числа элементов в системе. Сколько бы элементов в ней ни было - десятки, сотни или миллионы - для статистического расчета достаточно около 100/Q вычислений функции y(X).

К оглавлению

5.4. Пример статистического расчета

В качестве примера рассмотрим результаты расчета системы, изображенной на рис. 5.5, но с параметрами p1= p2= p3= 0.001, p4= p6= p9= 0.002, p5= p7= p8= 0.003, p10 = ... = p15 = 0.005.

В эксперименте, состоявшем из N = 200 000 опытов, отказное состояние системы осуществилось 35 раз. Отсюда Q` = 35/200000 =1.75 ·10-5.

Точность этой оценки в смысле коэффициента вариации можно найти с помощью (5.12)

(5.12)

или

Отсюда, используя правило трех сигм, получаем

Q = (1,75 ± 0,9) · 10-4

К оглавлению

5.5. Сравнение аналитического и статистического методов

Как аналитический расчет надежности, так и статистический имеют свои достоинства и недостатки.

Достоинство аналитического метода заключается в его аналитичности: он позволяет выразить искомую характеристику Q через заданные параметры p1,...,pn в общем виде (в виде формулы). Основной недостаток аналитического метода - трудоемкость расчета надежности неприводимых систем. Время расчета неприводимых систем становится пропорционально n · 2n (где n - число элементов системы), тогда как для приводимых систем - просто пропорционально n.

Основное достоинство статистического метода - как раз в возможности быстро рассчитать надежность любых систем, в том числе неприводимых. Трудоемкость расчета любых систем пропорциональна n.

Основной недостаток - статистический расчет получается приближенным и только для конкретных числовых значений p1 , ... ,pn. Если нужно исследовать зависимость Q от этих параметров, то приходится повторять расчет для их новых значений.

Значительное затруднение возникает и при расчете статистическим методом высоконадежных систем. Из (5.11) видно, что при малой вероятности Q число опытов N может оказаться слишком большим. Например, для системы на рис. 5.5. вероятность отказа Q имеет порядок 10-8 . Согласно (5.11) для расчета Q статистическим методом потребовалось бы порядка N =1010 опытов. Эксперименты с таким N трудно осуществлять даже на мощных современных компьютерах.

Таким образом, неприводимые высоконадежные системы - это особенно твердые орешки: они не рассчитываются ни аналитическим, ни статистическим методом.

В следующем разделе рассматривается аналитико-статистический подход, который соединяет достоинства как аналитического, так и статистического методов.

К оглавлению


6. АНАЛИТИКО-СТАТИСТИЧЕСКИЙ ПОДХОД

6.1. Общая характеристика
6.2. Метод взвешивания
6.3. Пример взвешенного моделирования
6.4. Метод расслоения
6.5. Примеры расслоения эксперимента

6.1. Общая характеристика

Под аналитико-статистическим подходом будем понимать изобретение всевозможных усовершенствований схемы статистического эксперимента с целью уменьшить число опытов N, необходимое для достижения заданной точности расчетов.

При этом усовершенствования схемы эксперимента выполняются путем ее формального, т.е. аналитического преобразования. Сам же эксперимент остается статистическим.

Здесь используется большое разнообразие идей; см., например, [5,9,10,12,16- 18,20,42].

Одна из наиболее известных - это идея понижения дисперсии. Она состоит в следующем.

Поскольку при заданной точности v число опытов N определяется формулой

то N можно уменьшить за счет снижения дисперсии D(y) выходной с.в. y. Точнее говоря, вместо y предлагается разыгрывать "немножко другую" с.в. y'с тем же значением искомого м.о. M(y') = M(y), но с меньшей дисперсией D(y') <<D(y).

Осуществление этой идеи рассмотрим на примере метода взвешивания и метода расслоения.

Система в целом также либо работает, либо находится в состоянии отказа. Состояние системы однозначно определяется совокупностью состояний элементов, т.е. является логической функцией вида:

К оглавлению

6.2. Метод взвешивания

Для простоты изложения метода рассмотрим статистический эксперимент с одномерной входной с.в. x (рис. 6.1).

Рис. 6.1. Исходная схема эксперимента

Формально искомое м.о. M(y) определяется выражением

M(y) = т y(t)f(t)Dt. (6.1)

Метод взвешивания основан на аналитическом преобразовании выражения (6.1) и на соответствующем изменении схемы эксперимента.

Разделим и умножим подынтегральное выражение в (6.1) на произвольную п.р.в. p(t), не равную нулю в пределах интегрирования. В результате мы увидим, что

(6.2)

где y'=y(x)f(x)/p(x), а x~p(t). Следовательно, вместо расчета оценки [image] можно выполнять расчет оценки [image]. Схема расчета [image] представлена на рис. 6.2. Такой переход от исходной схемы рис. 6.1 к схеме рис. 6.2 называется методом взвешивания. Обе схемы эквивалентны с точки зрения м.о. выходной с.в., но различны с точки зрения ее дисперсии.

Описанный одномерный вариант метода легко распространяется и на случай многомерной с.в. Для этого достаточно в (6.2) вместо скаляров x, t записать соответственно векторы X = (x1, ... , xn), T = (t1, ... , tn ). Обобщение на случай дискретной с.в. X осуществляется путем замены интеграла на сумму, а плотностей - на вероятности.

Искусство применения метода взвешивания сводится к подбору такой функции p, чтобы было D(y')<<D(y). Тогда получится N'<<N. Можно, например, подбирать п.р.в. p так, чтобы с.в. y'= y f / p была по возможности ближе к константе, хотя бы в каком-нибудь интуитивном смысле, так как при y'=сonst должно быть D(y') = 0.

Рис.6.2. Преобразованная схема

Основное достоинство метода взвешивания - в простоте преобразования схемы эксперимента. Программа вычисления функции y(x), т.е. модель исследуемой системы, не меняется. Добавляется только расчет множителя f/p (веса), да генератор входной с.в. с распределением f заменяется на генератор с распределением p.

К оглавлению

6.3. Пример взвешенного моделирования

Рассчитаем надежность системы, изображенной на рис. 5.5. Указанные для нее вероятности отказа элементов весьма малы, поэтому непосредственный статистический расчет надежности требует слишком больших затрат машинного времени (требуется порядка 1010 опытов).

Применяя в данной задаче метод взвешивания, получаем схему эксперимента, изображенную на рис. 6.3.

Рис. 6.3. Взвешенное моделирование надежности

Здесь P(X) определяется вероятностями отказа элементов p1,..., p15, приведенными на рис. 5.5. Вероятность P(X) конкретного значения с.в. X = (x1,..., x15) можно записать в виде

(6.3)

где qj =1 - pj. Запись (6.3) означает, что P(X) - это произведение, составленное из pj и qj, причем элементу xj=1 соответствует сомножитель pj, а xj= 0 - сомножитель qj. Например, набор X=(1,1,0,0,0,0,0,0,0,0,0,0,0,1,0) имеет вероятность P(X) = p1p2q3q4q5q6q7q8q9q10q11q12q13p14q15.

Распределение P’(X) может быть любым, лишь бы ненулевым P(X) соответствовали также ненулевые P’(X). Зададим P’(X) аналогично P(X), но изменив вероятности отказа элементов pj. Примем, что P’(X) задается вероятностями pўj = 4000pj. Таким образом, имеем pў1= pў2= pў3= 0.04, pў4= pў6= pў9= 0.08, pў5= pў7= pў8= 0.12, pў10= ... = pў15 = 0.2. Мы увеличили вероятность отказа системы, т.е. того, что y = 1. Однако благодаря весовому коэффициенту P(X)/P’(X) среднее значение выходной с.в. y’ остается равным Q.

При программировании эксперимента по схеме рис. 6.3 следует уделить достаточное внимание процедуре расчета веса P(X)/P’(X), которая может включать операции умножения на очень большие или на очень малые числа (для pўi умножений на такие числа не возникает).

Описанный выше конкретный вариант взвешивания позволяет оценить вероятность отказа системы электрообеспечения при не слишком больших N. Так, для N = 1000 в эксперименте были получены следующие оценки:

, [image].

Таким образом, используя правило трех сигм, для исследуемой системы можно записать

Q = (1,34 + 0,66) · 10-8. (6.4)

К оглавлению

6.4. Метод расслоения

Для типовой исходной схемы статистического эксперимента (рис.6.4) введем следующие обозначения:

W - пространство исходов, т.е. множество всех возможных значений с.в. X.

W1, ...,W k слои, т.е. такие подмножества пространства W , что

wj = P {XОWj} - вероятность j-го слоя.

Mj - условное математическое ожидание с.в. y в слое x О W j.

Рис.6.4. Исходная схема эксперимента

Известно, что безусловное м.о. с.в. выражается через условные следующим образом:

(6.5)

Отсюда происходит идея расслоения эксперимента: вначале выполнить эксперименты в отдельных слоях X О W j, в которых с.в. y более или менее "близка к константе", а затем по найденным "почти точным" оценкам [image] вычислить в соответствии с (6.5) оценку

(6.6)

Схема расслоенного статистического эксперимента, воплощающая эту идею, представлена на рис. 6.5.

Для контроля за точностью эксперимента можно в каждом слое вычислять условную дисперсию выходной с.в.:

X О W j (6.7)

и дисперсию внутрислойной оценки

(6.8)

Эксперимент в каждом слое проводится точно также, как по исходной схеме. Только вместо безусловной п.р.в. f(t) используется условная fj(t):

Рис.6.5. Схема расслоенного эксперимента

В соответствии с линейным выражением (6.6), дисперсия итоговой оценки может быть вычислена по формуле:

(6.9)

Если нужно знать ошибку итоговой оценки в форме коэффициента вариации, то она вычисляется как обычно:

(6.10)

Следует иметь в виду, что кроме удачного выбора слоев существует дополнительная возможность выигрыша в точности - за счет правильного распределения общего числа опытов N по слоям: N = N1 + ... + Nk. Опыты можно распределять рационально - в соответствии с формулой

j =1,…,k, (6.11)

или оптимально - по формуле

j =1,…,k, (6.12)

где [image] - условное среднеквадратическое отклонение с.в. y в слое W j. Обычно sj априорно неизвестны и их заменяют на оценки

, j =1,…,k, (6.13)

где [image] находят по (6.7) путем выполнения пробного прогона программы с не очень большими Nj.

При использовании распределения опытов (6.12) ошибка v расслоенного эксперимента будет минимальной для данного N. Распределение (6.11) приводит к худшей (во всяком случае не лучшей) точности, чем (6.12). Тем не менее, распределение опытов (6.11) используется чаще, так как оно проще и все же гарантирует (при любом расслоении эксперимента), что ошибка получится меньше (в крайнем случае - не больше), чем в исходной схеме. Эти свойства метода расслоения доказываются достаточно просто.

Таким образом рациональное, а тем более и оптимальное распределение опытов гарантирует, что даже при самом неудачном расслоении мы не ухудшим точности расчетов.

К оглавлению

6.5. Примеры расслоения эксперимента

Первый пример

Рассмотрим методику применения метода расслоения к расчету высоконадежных структурно сложных систем на примере системы, изображенной на рис 5.5.

Пространство W исходов в данном случае содержит 215 значений двоичной с.в. X =(x1,...,x15). Обозначим через s число отказавших элементов, т.е. s= x1,..., x15= s(X) и разобьем W на слои W0, ... ,W3 следующим образом:

W0= {X : s (X) = 0},

W1= {X : s (X) = 1},

W2= {X : s (X) = 2},

W3= {X : s (X) і 3},
(6.14)

Выразим безусловное м.о. M(y) через условные м.о. в слоях W0, ... ,W3:

(6.15)

где wj = P {Wj}, j=0, ... ,3. Заметим, что, как видно из рис. 5.5, при числе отказавших элементов 0 или 1 система не отказывает. Следовательно, в слоях W0 и W1 всегда y=0 и, значит, M0= M1= 0. Таким образом, равенство (6.15) принимает вид

M(y) = w2M2 + w3M3 (6.16)

Кроме того, вследствие малости w3 слагаемое w3M3 также можно отбросить. Тогда из (6.16) с точностью до 2-3-х значащих цифр имеем

M(y) = w2M2 (6.17)

где w2 = P{W2} = 1.0976 · 10-7и весь эксперимент можно проводить только в слое W2.

Схема расслоенного таким образом эксперимента приведена на рис.6.6.

Рис.6.6. Схема расслоенного эксперимента для системы (рис.6.5)

На этой схеме P2 - условное распределение с.в. X = (x1, ... x15) для условия s(X)=2. Генерация с.в. X~P2 может осуществляться разными методами. Они рассматриваются ниже.

По схеме рис. 6.6 был выполнен эксперимент с числом опытов N=1000. В результате в 150 случаях был обнаружен отказ системы. Следовательно

= 150/1000 = 0.15, (6.18)

(6.19)

Определим дисперсию оценки (6.18) с учетом правила трех сигм

M(y) = (1.65 ± 0.35) 10-8. (6.20)

Второй пример

Попробуем еще один вариант расслоения: составим каждый слой ровно из одного исхода с двумя или тремя отказами элементов. Тогда мы имеем всего 560 таких слоев:

      W 1= {X: X = (1,1,0,0,0,0,0,0,0,0,0,0,0,0,0)},
      .
      .
      .
      W 560 = {X: X = (0,0,0,0,0,0,0,0,0,0,0,0,1,1,1)}.

Ясно, что в каждом таком слое имеет смысл производить только по одному опыту, и общее число опытов получается небольшим: N = 560. С.в. X в каждом слое вырождается в одно детерминированное значение.

В результате прогона модели с таким расслоением была получена оценка:

(6.21)

Все приведенные здесь цифры - точные. В целях контроля выполнен также расчет системы при p1=... =p15= 10-6. Для таких одинаковых значений pi известно точное аналитическое решение

M(y) = 1.12 · 10-8. (6.22)

Статистический эксперимент с 560 слоями и одним опытом в каждом слое дал оценку

. (6.23)

Как разыграть состояния элементов?

Генерация с.в. X ~ P2 должна выполняться в соответствии с условным распределением вероятностей

P2 = P {X/s (X)=2}. (6.24)

Здесь ненулевую вероятность имеют только значения X с двумя единичками. Если эти единички находятся на местах i, j (что соответствует отказу i-го и j-го элементов, i , j = 1...n), то

(6.25)

Разделив последнее выражение на q1...qn, получим более простую запись:

l, m О{1,..., n} (6.26)

Теперь введем обозначения:

(6.27)

Тогда окончательно получаем

(6.28)

где i , j - номера единичных элементов вектора X О W2.

Один из методов реализации X состоит в разыгрывании пары (i, j) ~ P(i, j) = rirj / S2. Для этого применяется схема разыгрывания случайных событий (см. п. 2.2). Этот метод, однако, удобен только при не очень больших n, так как иначе число возможных событий (число пар i, j), равное n(n+1)/2, может оказаться слишком большим.

Другой способ реализации с.в. X состоит в том, чтобы построить процедуру, состоящую в разыгрывании сначала номера i одного элемента, потом номера j второго элемента. При этом результирующая вероятность выбора пары (i, j) должна совпадать с заданной в (6.28). Построение такой процедуры - задача не тривиальная. Одна из наиболее простых процедур такого типа представляется в виде следующего алгоритма, написанного с учетом обозначений (6.27).

(1) Реализуем с.в. i ~P (i) = 0.5 ri (S1- ri) / S2; i=1,...,n.

(2) Реализуем с.в. j ~P (j/i)=rj / (S1- ri); j = 1,...,n, j№i

(здесь i - значение, выбранное на шаге (1)).

Докажем, что вероятность P(i, j) выбора пары i, j для этого алгоритма соответствует требуемому закону (6.28): -

P(i,j)=P(i)P(j/i)+P(j)P(i/j) = [0.5ri(S1- ri)/S2][rj/(S1- ri)]+[0.5rj(S1- rj)/S2][ri/(S1- rj)]=
=0.5 ri rj/S2+ 0.5 rj ri /S2= ri ri /S2, ч. и т.д.

Также нетрудно убедиться и в том, что на обоих шагах алгоритма мы имеем дело действительно с распределениями вероятностей, т.е. что P(i) и P(j/i), j № i - неотрицательные величины, дающие в сумме единицу:

.

Предельно простой и универсальный метод генерации с.в. X внутри любых слоев типа S(X)=1, S(X)=2 и т.д. изложен, пожалуй, Б.И.Плаксом в [20]. Идея этого метода состоит в том, чтобы вначале перейти к взвешенной схеме моделирования (см. пп. 6.2 и 6.3), в которой вероятности отказа всех элементов pi сделать одинаковыми, а затем уже взвешенную систему моделировать путем расслоения эксперимента. Тогда в силу равенства всех pi любые номера отказавших элементов в слое равновероятны и их выбор может осуществляться просто по (2.2). Так, в слое S(X)=2 для системы рис. 5.5 вначале будет вычисляться i по формуле

i = л 15 z ы + 1,

где z - значение БСВ, а затем j по такой же формуле. Если j совпадет с i, то вычисление j следует повторить.

Кроме того, в методе Плакса очень легко рассчитываются вероятности слоев, используемые в (6.15).

Б.И.Плакс также оценивал своим методом надежность системы, воспроизведенной на рис. 5.5, при N = 1000 и нашел, что

Эффективность метода Плакса детально исследована в [18]. В ней же путем усовершенствования метода и был получен результат (6.21).

К оглавлению


7. РЕАЛИЗАЦИИ АНАЛИТИКО-СТАТИСТИЧЕСКОГО ПОДХОДА

Назад

Содержание

Вперед

7. РЕАЛИЗАЦИИ АНАЛИТИКО-СТАТИСТИЧЕСКОГО ПОДХОДА ПРИМЕНИТЕЛЬНО К АНАЛИЗУ МОДЕЛЕЙ ТЕЛЕКОММУНИКАЦИОННЫХ СЕТЕЙ

7.1. Задачи системного проектирования ТКС
7.2. Статистическое оценивание функциональной надежности сети
7. 3. Аналитико-статистический расчет вероятностно-временных характеристик транспортировки пакетов по сети
7.4. Модель и методика расчета малых вероятностей потерь в буфере конечной емкости

7.1. Задачи системного проектирования ТКС

Этап системного проектирования ТКС включает структуризацию сети, построение сетевой модели, построение функциональных моделей элементов сети, анализ структурно-функциональной модели сети, заключающийся в определении вероятностно-временных характеристик (ВВХ) сети в зависимости от параметров отдельных функциональных моделей и организации связей между ними.

Рассматривая ТКС как средство взаимодействия множества пользователей (включая и автоматические, например, ЭВМ) между собой, можно выделить три группы задач, решаемых на этапе системного проектирования.

Первая группа задач связана с синтезом топологической структуры ТКС с учетом ее развития, динамики потоков и живучести, выбором производительности и алгоритмов функционирования ее элементов (узлов коммутации и каналов связи).

Вторая группа задач обусловлена реализацией технологии взаимодействия, которая получила свое оформление в виде системы протоколов обмена информацией, имеющих многоуровневую структуру.

Третья группа задач - это задачи управления взаимодействием на уровне транспортной сети; решение их связано с разработкой эффективных процедур управления коммутацией, маршрутизацией, управлением потоками, реконфигурацией сети.

При системном проектировании сети в качестве критериев используются ВВХ доставки сообщения, стоимость сети, структурная надежность либо живучесть сети. Причем стоимость сети выступает часто в виде ограничения. Названные задачи системного проектирования связаны между собой. Желательно комплексное решение, что является сложной задачей. На практике сложился подход, когда решаются отдельные задачи и затем объединяются в общей процедуре проектирования. На первом этапе эвристическим методом определяется исходный граф сети. Затем проводится анализ исходного варианта с привлечением аналитических и имитационного методов моделирования. Определяются надежностные и временные характеристики ТКС. Определяется область допустимых значений нагрузки, при которых обеспечивается требуемое время передачи информации. Сложность задачи заключается в том, что для различных участков сети должны выполняться различные как в качественном, так и в количественном отношении условия транспортировки пакетов. Поэтому в большинстве случаев в роли критериев эффективности может выступать время доставки пакетов между отдельными корреспондирующими парами "источник-адресат", что обусловливает необходимость оценки ВВХ доставки между всеми корреспондирующими парами узлов при заданных структуре сети, алгоритмах маршрутизации и дисциплине очередей.

В этом разделе мы продемонстрируем использование аналитико-статистического метода и построение соответствующих моделей применительно к решению трех задач: оцениванию функциональной надежности сети, оцениванию времен доставки пакетов по сети и расчету необходимой емкости буферных накопителей, применяемых в узлах коммутации.

Внимание к расчету емкости буферного накопителя обусловлено особенностями сетей АТМ.

В сетях АТМ управление потоком не применяется, так как его трудно осуществить при высоких скоростях передачи. Обслуживание в таких сетях характеризуется распределением ресурсов. Статистическое распределение ресурсов определяется с учетом вида трафика с тем, чтобы обеспечить низкие вероятности потерь пакетов данных и речи и незначительный разброс задержек пакетов при передаче видеоинформации.

К оглавлению

7.2. Статистическое оценивание функциональной надежности сети

Характеристики функциональной надежности отражают возможность сети предоставлять связь абонентам в течение заданного интервала времени начиная с произвольного момента - момента поступления вызова (требования). Другими словами, функциональная надежность сети во многом определяется процессом установления соединения на сети. Установление требуемого соединения на сети, несмотря на повреждения или отказы отдельных элементов, позволяет сохранять работоспособность, т.е. обеспечивает живучесть сети. Для оценки функциональной надежности сети количественным показателем может служить вероятность установления (неустановления) соединения на сети при поступлении соответствующего требования. В свою очередь процесс установления соединения зависит от принятой системы назначения маршрутов и распределения потоков информации.

Время установления соединения жестко ограничивается, как и вероятность превышения этого времени; ограничивается также вероятность неустановления соединения за допустимое время.

Практически любая стратегия маршрутизации представляет собой совокупность маршрутных таблиц (МТ), которые распределяются по узлам коммутации и указывают, как в зависимости от конечного адресата должен быть распределен по выходным линиям трафик, поступающий в данный узел. Такая стратегия обеспечивает ряд альтернативных путей (направление обмена) для каждой пары корреспондирующих абонентов. Получаемое статическое решение МТ подлежит оценке с учетом динамики нагрузки и проверяется средствами имитационного моделирования. В качестве такого средства можно использовать имитационную модель установления соединения. Необходимость названной проверки обусловлена ответственностью получаемого решения - МТ обеспечивают живучесть сети.

7.2.1. Концепция моделирования установления соединения

Процесс установления соединения представляет собой посылку и прохождение вызова от источника к адресату по путям, составляющим направление обмена между источником и адресатом. Направление обмена есть множество путей - "гамак" кратчайших путей (ГКП), соединяющий источник вызова с адресатом (рис. 7.1,а,б).

В общем случае сеть неоднородна, она может содержать два типа станций (узлов): станции, допускающие альтернативные исходящие пучки каналов, и станции, которые допускают только один исходящий пучок каналов в соответствии с номером, содержащимся в поле вызова.

Пути между станциями состоят из цепочки транзитов. В свою очередь каждый транзит (ветвь) представляет собой пучок каналов, соединяющих два смежных узла в такой цепочке. Взаимное сочетание названных типов станций обусловливает особенности и влияет на характеристики процесса доставки вызова: если в транзите все каналы требуемой скорости заняты, то вызов через такой транзит не пойдет и должен быть передан по другому пути, т.е. за счет рестарта возвращается либо к одной из предшествующих станций первого типа, либо к станции-источнику для повторной посылки. Эта особенность должна быть отражена в модели.

Кроме фазы доставки вызова до адресата процесс установления сквозного соединения в режиме "обратной волны" включает фазу проключения коммутационных полей станций, через которые вызов прошел при его доставке к конечной станции, и тем самым в промежуточных станциях был зарегистрирован необходимый для сквозного соединения ресурс.

Рис. 7. 1. Направление (а) и гамак кратчайших путей (б) между истоком i и стоком j

Таким образом, с точки зрения прохождения вызова от станции к станции состояние транзита (ветви) бинарно: либо транзит полностью занят и вызов по нему не проходит (состояние транзита "1"), либо транзит имеет хотя бы один незанятый канал требуемой скорости, и вызов через такой транзит пройдет (состояние транзита "0"). Соответственно и конечный результат установления соединения имеет два исхода: либо вызов дойдет до адресата и в "обратной волне" проключением коммутационных полей станций будет установлено соединение за суммарное время (доставка вызова и проключение полей), не превышающее допустимого, либо соединение не будет установлено, поскольку вызов вообще не дойдет до адресата (во всех путях "гамака" наличествуют полностью занятые нагрузкой транзиты), либо дойдет до адресата, но суммарное время превысит за счет рестартов и повторных попыток доставки вызова по другим путям.

Проключение коммутационных полей в режиме "обратной волны" - детерминированный процесс: если вызов доставлен до адресата, то в рамках поставленной задачи проключение происходит всегда и практически за фиксированное время. Фактор случайности связан с фазой доставки вызова.

Факт доставки вызова устанавливается с помощью имитации процесса доставки. При этом ввиду возможности рестартов - возвращения вызова на предыдущие узлы - число транзитов, пройденных вызовом при его доставке к адресату, оказывается случайным. Именно этот аспект придает стохастичность процессу установления соединения. Времена обработки вызова при его прохождении через транзит, времена проключения коммутационных полей в режиме "обратной волны" и, наконец, времена рестартов можно считать фиксированными величинами. Тогда время доставки вызова - случайная величина - определяется выражением

tУС = ( ngtтр + nрtов + rtr) (7.1)

где ng- число транзитов, пройденных вызовом при его доставке адресату с учетом случившихся рестартов; nр- число транзитов в пути, по которому устанавливается соединение в режиме "обратной волны"; r - число рестартов, имевших место в зафиксированной реализации процесса установления соединения. В общем случае 0 Ј r < rДОП ,( rДОП - допустимое число рестартов при установлении соединения). При r >rДОП считается, что вызов адресату не доставлен; tтр,tов,tr - времена прохождения транзита, проключения коммутационного поля в узле, рестарта, соответственно.

При имитации доставки вызова каждый прогон модели в случае доставки вызова до адресата дает реализации случайных величин ng, nри r, по значениям которых в соответствии с выражением (7.1) подсчитывается время доставки вызова tУС. Если tУС Ј tДОП, то в результате испытания (одного прогона модели) имеем факт установления соединения. Если же при доставке вызова tУС > tДОП, либо при имитации доставки оказалось r > r ДОП, либо вообще в разыгранной реализации "гамака" не оказалось путей со свободными каналами требуемой скорости, то в результате испытания имеем факт неустановления соединения.

Входными данными для задачи являются структурная и потоковая метрики, гамаки кратчайших путей, полученные на их основе значения вероятностей полной занятости транзитов (ветвей) и ограничения, при которых должно выполняться соединение: максимально допустимое число транзитов в соединении, допустимое время установления соединения (либо доставки вызова), времена продвижения вызова через транзит и проключения коммутационного поля, число рестартов (возвратов вызова для передачи его по обходным путям, время выполнения рестарта).

Выходными данными являются статистики, позволяющие оценить вероятность установления соединения (доставки вызова адресату) за время, не превышающее допустимого, и средние и среднеквадратические значения tУС, характеризующие процесс установления соединения, которые могут быть использованы при упрощенных аналитических расчетах.

7.2.2. Модель процесса установления соединения

С точки зрения установления соединения направление обмена пары станций (i, j) можно представить в виде графа (рис.7.1), G = (A,L), где A = (F, S) - множество узлов (станций), входящих в направление обмена; F = {Fk}, k = [image] - множество узлов, отображающих станции 1-го типа, допускающие альтернативные исходящие пути; S = {Sw}, [image] - множество узлов, отображающих станции 2-го типа, допускающие единственный исходящий путь в соответствии с адресом вызова; L={hn}, [image] - множество ветвей, отображающих линии связи (транзиты), соединяющие узлы направления обмена; hn = (n1,...,nc) - упорядоченный вектор скоростей, на которых может осуществляться связь по данной hn ветви (транзиту). Упорядоченность заключается в том, что n m' > n m при m < m', m, m'= [image] .

Полагаем, что установление соединения с требуемой скоростью n m может быть осуществлено, если в пучках транзитов пути, по которому доставлен вызов до адресата, были свободными каналы со скоростями n m' і n m, m = [image], m' = [image] .

В качестве критерия эффективности ТМ примем вероятность установления соединения за время, не превышающее допустимого. В общем случае, в силу большой размерности сети для решения задачи оценки вероятности установления (неустановления) соединения на сети целесообразно использовать двухступенчатую выборку.

На первой ступени из генеральной совокупности станций сети {Ai}, i = [image] производится случайная выборка объемом N 1 пар станций (i, j), i, j = [image], i № j , для которой должна быть оценена возможность установления соединения. Пара (i, j), в которой i - источник вызова, а j - адресат, задает для второй ступени входные данные, которыми являются записи файлов, содержащие состав ГКП i j и соответствующие вероятности, характеризующие "закрытость" транзитов гамака.

При случайном выборе пары (i, j) в начале в соответствии с вероятностями активности станций p(i), i = [image] выбирается номер станции-источник вызова, а затем в соответствии с условными вероятностями активности связей источника p(j Ѕ i), j = [image], i № j разыгрывается номер станции-адресата.

Файлы вероятностей p(i) и p(j Ѕ i), i,j = [image], i № j организуются при решении задачи распределения потоков, причем

и [image] ,

где lik - интенсивность нагрузки, поступающей от k-го абонента, подключенного к i-й станции-источнику; lij - интенсивность обмена информацией между источником i и адресатом j, задаваемая матрицей нагрузки.

Выходом цикла являются номера пар станций (i, j), где i - источник, j - адресат, и значение вероятности (если требуется) p(i,j) = p(i)p(j Ѕ i).

Отметим, что для ускорения моделирования на первой выборочной ступени все множество гамаков {Гij}, i = [image], i№ j можно декомпозировать на подмножества, используя подход, изложенный в р. 7.3. Для объединения гамаков в подмножества достаточно учитывать близость условий передачи в основных путях гамаков.

На второй ступени производится расслоение по числу h транзитов (h = hmin... hmax) в Гij, полностью занятых нагрузкой при требуемой скорости. Случайным образом в отдельном опыте разыгрывают h номеров полностью занятых транзитов в Гij), и на полученной реализации Гijимитируется процесс доставки вызова адресату. По факту доставки вызова фиксируются значения ng, np и r . Процедура при случайном выборе h полностью занятых транзитов повторяется N 2 раз. По результатам экспериментов вычисляются необходимые статистики.

Благодаря равновероятности генерирование реализации Xh предельно упрощается и сводится к "выбору наугад" номеров ветвей (транзитов) из общего числа d, закрытых для вызова. Очередной номер определяется по формуле j = [Ud] + 1, где U - случайное число, U О[0,1], получаемое путем обращения к датчику случайных чисел. Скобки [ ] означают округление в меньшую сторону. Составляющие гамака an с этим номером полагаются равными единице.

7.2.3. Аналитико-статистический расчет оценки вероятности
установления (неустановления) соединения на сети

Так же, как и при оценке вероятности связности сети, задача расчета характеристик установления- неустановления соединения может быть сведена к задаче оценивания математического ожидания M x СВx = f(an), причем СВ an=(a1,... , ad) имеет закон распределения вероятностей p (т.е. p), который известен. В данной задаче xО{0,1}, x = 0 соответствует установлению соединения, =1 - неустановлению соединения. M x = P{x=1} имеет смысл вероятности неустановления соединения. СВx невырождена, т.е. 0 < M x < 1 .

СВan интерпретируем как вектор an = ( [image] ), отображающий состояние ветвей (транзитов), входящих в моделируемый гамак, с точки зрения наличия или отсутствия свободных каналов со скоростью u , которая требуется (допустима) для устанавливаемого соединения. СВ [image] О{0,1}, i = [image], d - число ветвей (транзитов) в моделируемом гамаке. Значение [image] = 1 понимается как отсутствие свободных каналов для скорости u в i-й ветви. СВ [image] считаются независимыми. В этом случае для каждого значения скорости u СВ a имеет конечное множество значений u О X, X = {xj; j = 1,... , n; n = 2d}. Для нее распределение (х), х О Х задается набором вероятностей p [image] . С учетом расслоения по h оценку Мx находим в виде:

(7.2)

где [image] - оценка вероятности неустановления соединения при h непроводящих вызов транзитах.

Оценка [image] находится методом РВМ. При выборе нового распределения считаем все возможные реализации гамака, содержащие h непроводящих ветвей (транзитов), равновероятными, т.е. примем вероятность реализации хО(X h) в виде

где w (x) - вес реализации хО(X h). Условие нормировки S q(x)=1 при этом выполняется. Новая осредняемая функция имеет вид

и оценка [image] x(xIh)k [image] ,

где [image] - число разыгранных состояний (реализаций) гамака Гij; x(xIh)k- исход доставки-недоставки вызова при k-м испытании (k-й реализации Гij); x(xIh)kО(0,1) .

Для задания состояния гамака Гij(h ) на k-м розыгрыше применяем случайный выбор h номеров непроводящих ветвей (транзитов). В результате вектор получает конкретную реализацию x, содержащую h единиц и d-h нулей. В соответствии с правилом прохождения вызова по гамаку путей Гij и правилом установления соединения вычисляем значение

и вероятность полученной реализации

.

Значение x(xIh)k= 1 имеет место, если вызов не дошел до адресата (все пути гамака оказались непроводящими), либо число рестартов превысило допустимое значение, либо вызов дошел до адресата и установлено соединение, но за время tУС > tДОП .

Таким образом, оценку вероятности неустановления соединения между парой станций (i, j) при h непроводящих ветвях (транзитах) в гамаке Гij определяем в виде

.

Вероятность p( [image] ) рассчитывается по формуле

.

Вероятности р( [image] =1), i = [image] заданы при решении задачи распределения потоков.

Окончательная оценка вероятности неустановления соединения между парой станций (i,j) имеет вид

Таким образом, при решении рассмотренной задачи применен метод расслоения и равномерно взвешенный метод генерации с.в. внутри слоев, предложенный Б.И.Плаксом.

К оглавлению

7. 3. Аналитико-статистический расчет
вероятностно-временных характеристик транспортировки пакетов по сети

При построении полной (топологически подобной) имитационной модели сети в целях оценки ее ВВХ возникает задача снижения размерности модели.

Формально уменьшение размерности может быть сведено к выделению на структуре сети некоторого множества классов объектов с последующей их заменой соответствующими моделями. Актуально при этом обеспечить однотипность моделей для представителей различных кластеров. Важно также, чтобы эти модели с достаточной полнотой отражали специфику процессов, определяющих ВВХ сети.

В качестве таких представителей различных кластеров целесообразно выбрать путь обмена информацией (ПОИ). Сеть декомпозируется на подмножества ПОИ. В каждое подмножество входят пути, близкие по условиям транспортировки пакетов. Моделью ПОИ и объектом имитации может служить виртуальный канал (ВК), представляющий собой некоторый маршрут в сети, состоящий из последовательности узлов коммутации (УК) и каналов связи (КС) от узла-источника до узла-адресата. Одна и та же модель ВК для представления разных кластеров будет отличаться только параметрами.

Для имитации взаимодействия с другими направлениями обмена информацией, т.е. учета влияния потоков, циркулирующих по сети и являющихся транзитными в выделенном ВК, строится вероятностный эквивалент нерассматриваемой части сети в виде генератора транзитных (фоновых) потоков.

Такой подход даст выигрыш в том смысле, что имитируемый объект будет существенно меньшей размерности. Как теоретическую основу данного подхода можно рассматривать метод расслоенной выборки, предлагая в качестве совокупности факторных признаков все множество ПОИ сети, поскольку передача информации по отдельным направлениям обмена или маршрутам происходит с некоторым дискретным распределением вероятностей использования маршрутов, а также функциями распределения поступления сообщений в сеть и распределениями длительностей обслуживания на узлах коммутации и каналах связи.

7.3.1. Декомпозиция сети на подмножестве путей обмена информации

Специфика декомпозиции ТКС, как мы отметили, заключается в том, чтобы при моделировании выделенной части воспроизвести взаимодействие процессов, протекающих во всей сети и влияющих на характеристики выделенной части. Это требование реализуется при декомпозиции ТКС на классы направлений обмена информацией. Кластеризацию направлений можно выполнить по признакам ПОИ, по которому проложен основной маршрут. При этом факторными признаками выступают длина пути (число транзитных участков в ПОИ) и нагрузка, приходящаяся на УК и КС, составляющие рассматриваемый путь. Поскольку значения нагрузки УК и КС путей могут в сильной степени отличаться друг от друга, необходимо проранжировать УК и КС сети: канал или узел будет обладать большим рангом, если на него приходится большая нагрузка. Определение значений нагрузки для каждого элемента сети можно осуществить в программе построения маршрутов, задав, кроме графа сети, производительности l i каналов и узлов. Тогда в соответствии с матрицей нагрузки для каждого элемента сети будут подсчитываться суммарная интенсивность l i входящего потока и значения нагрузки r i = l i/m i. На полученном множестве значений нагрузки узлов и каналов каждому элементу назначается ранг. Подробнее на вопросе назначения рангов остановимся ниже.

Итак, пусть m1- количество рангов для УК, m2- количество рангов для КС. Тогда для каждого пути прохождения информации длиной в d транзитов набор признаков есть вектор Z = (z(1),...,z(2d+1)), компоненты которого [image] принимают определенные значения в зависимости от рангов элементов, составляющих путь. Размерность вектора признаков 2 d+1 (d+1 - число УК, d - число КС) меняется в зависимости от длины пути, т.е. числа транзитных участков.

Разбиение множества путей П={p1,...,pNd} проводится с помощью двухэтапной классификации объектов. На первом этапе по признаку, характеризующему длину пути, образуем подмножества в соответствии с градациями этого признака. Получим D подмножеств Пd, [image] , где D - максимальное число транзитных участков для путей piО П, [image]. Далее, внутри каждого подмножества Пd проводим разбиение на классы по степени значимости, используя процедуру кластерного анализа. В итоге любой объект piО П(d) характеризуется набором признаков zi = (zj(1),..., zj(2d+1)), где компоненты [image] , [image] имеют диапазон изменения z1 Ј [image] Ј [image], [image] - значение максимального ранга, назначенного УК; компоненты [image] , [image] имеют диапазон изменения z1 Ј [image] Ј [image] , [image] - значение максимального ранга, назначенного КС. При этом элемент сети при наличии большей нагрузки обладает большим значением ранга.

Группирование путей pjО П(d)в непустые непересекающиеся подмножества, именуемые классами, осуществляется по критерию близости объектов. Мера близости представляет собой функцию, ставящую в соответствие каждой паре точек некоторое число Sij, характеризующее степень сходства (близости) между объектами piи pj. Для нашего случая в качестве меры близости целесообразно принять функцию

(7.3)

Здесь S'ij- расстояние между точками pi и pj, определяемое в виде

где b - величина, характеризующая степень различия элементов, обладающих одинаковым рангом; величина (Zi(k) - Zj(k) )b выражает различие по k-му признаку объектов pi и pj.

С помощью меры близости (7.3) в один класс объединяются схожие между собой объекты, причем степень сходства у объектов, принадлежащих к одному классу, должна быть больше, чем степень сходства между объектами, относящимися к разным классам.

Для построения оптимального разбиения множества ПОИ в качестве критерия оптимальности классификации используем сумму "внутренних" связей за вычетом некоторого порогового значения, характеризующего существенность связей. Оптимальная по заданному критерию классификация удовлетворяет требованию - средняя связь внутри каждого класса превышает среднюю связь объектов этого класса как со всеми остальными классами объектов, так и с отдельными объектами каждого другого класса.

Пусть а - число, выражающее уровень существенности связи так, что связь Sij > а является существенной, и (Sij - а)> 0 указывает на сходство объектов, тогда как связь Sij< 0 является несущественной и величина | Sij - а| при Sij< 0 характеризует степень различия между объектами pi и pj

Классификация объектов характеризуется разбиением R = {R1,..,RL} множества объектов на непересекающиеся классы R1,..,RL. Показателем качества разбиения R = {R1,..,RL} при значении а порога существенности является число:

(7.4)

Иначе говоря, в качестве показателя связи между pi и pj принимается (Sij- а); связи внутри класса Rh характеризуются суммой [image] так, что f(а, R) есть сумма связей (Sij - а) внутри разбиения R. Оптимальным будет то разбиение, которое максимизирует f(а,R) при заданном а по всем возможным разбиениям данного множества объектов.

Матрица S*, определяемая мерой близости (7.3), симметрична, и ее элементы принимают значения в диапазоне от 0 до 1, т.е. Sij = Sji, 0Ј Sij Ј 1. Условимся, что связи Sii объектов с самими собой не рассматриваются, так что f({pi },{pi })=0 по определению.

Критерий оптимальности (7.4) позволяет проводить разбиение на априорно не заданное число классов. К достоинствам показателя (7.4) относится также возможность выбирать значение порога существенности связей на основе содержательного смысла конкретной задачи. При различных значениях порога оптимальными являются, вообще говоря, разные разбиения. В задаче классификации путей прохождения информации от значения порога зависит, в конечном счете, точность оценки для сети в целом.

7.3.2. Назначение рангов УК и КС.
Выбор порога существенности связей

Очевидно, что элементы сети, на которые приходится в общем-то разная нагрузка, могут иметь одно и тоже значение на введенной шкале расстояний между объектами, т.е. обладать одним и тем же значением ранга. Проблема в том, как определить диапазон изменения нагрузки, чтобы элементам сети, с разностью между их значениями нагрузок, не выходящими за пределы этого диапазона, можно было назначить один и тот же ранг.

При решении этой проблемы исходим из цели задачи. Цель состоит в проведении классификации объектов множества П= p1,...,pNd по набору факторных признаков таким образом, чтобы результативные признаки объектов, входящих в один класс, отличались между собой меньше, чем результативные признаки объектов, принадлежащих разным классам.

Нагрузка, приходящаяся на элементы сети, является факторным признаком. Точная зависимость исследуемой характеристики от факторных признаков неизвестна, иначе не было бы необходимости в имитационном моделировании. Но можно аппроксимировать эту зависимость приближенно, учитывая тот факт, что каждый элемент объекта (в нашем случае УК, КС) вносит свое влияние на исследуемую характеристику. Эту "долю влияния" назовем составляющей результативного признака (СРП). В свою очередь, по построенной приближенной зависимости искомой характеристики от факторных признаков можно назначить ранги элементам сети, а затем определить и значение порога существенности а.

Пусть Ti = j (ri) выражает зависимость СРП от фактора ri(нагрузки элемента i).Тогда, задав число b, как допустимое различие значений СРП в процентах либо долях, можно элементам i и j присвоить одинаковый ранг при выполнении условия

Ѕ1 - Ti/TjЅЈ b (7.5)

Рассматривая УК (КС) как систему MЅ DЅ 1 - на входе пуассоновский поток с интенсивностью l сообщений в единицу времени и детерминированное время обслуживания, можно оценить СРП в виде:

где [image] - постоянное время обслуживания.

Ранжирование элементов сети производится раздельно, сначала для УК, затем для КС по следующей схеме:

  1. Выбирается УК (КС), имеющий наименьшую нагрузку, проверяется выполнение условия (7.5) при сопоставлении характеристики выбранного узла (канала) с характеристиками всех остальных. Всем УК (КС), для которых это условие выполняется, присваивается ранг 1.

  2. Из оставшихся УК (КС) выбирается узел (канал) с наименьшей нагрузкой [image] (q=2 ) и выполняется процедура первого этапа. Группе выбранных УК (КС) присваивается ранг

где й x щ - ближайшее целое, не меньшее х.

    3. Процедура повторяется для q = 3,4,... до тех пор, пока не будет исчерпано все множество УК (КС).

Таким образом, ранжирование элементов определит диапазон изменения признаков Z(k); k = 1,...,2d + 1. Число рангов m зависит от задания величин значения b. При малом b число рангов больше, и наоборот, при увеличении b число рангов уменьшается. Если Zmax- максимальное значение ранга, назначенное УК (КС), а m - количество рангов, которое присвоено УК (КС), то в общем случае Zmaxі m .

Объект [image] О П(d)имеет 2 d+1 элементов (d-КС и (d+1)-УК), каждый из которых характеризуется некоторым рангом. Полное совпадение рангов элементов двух объектов [image] и [image] выражается совпадением значений векторов признаков Zi и Zj и в соответствии с (7.2) значение меры близости Sij = 1. В этом случае разница e результативных признаков объектов [image] и [image] будет лежать в пределах 0... (2 d+1)b. При полном расхождении векторов признаков двух объектов разница e результативных признаков не превышает величины

[Zm1(d + 1) + Zm2 d]b

где [image] , [image] - максимальные значения рангов УК и КС соответственно.

Значения функции Sij изменяются в пределах от 0 до 1. Чем больше сходство объектов, тем ближе значение к единице; чем больше разница результативных признаков, тем ближе значение Sij к нулю.

Принимая во внимание, что область изменения функции Sij от разницы результативных признаков невелика, и наблюдается обратная зависимость Sij от разницы результативных признаков, можно аппроксимировать изменение функции Sij (e) = S(e) линейной зависимостью (рис.7.2)

Задавая eдоп- желательную степень однородности объектов внутри классов, по (7.6) можно найти значение меры близости eдоп S(eдоп) = а, при котором достигается заданная степень однородности, т.е. расхождение результативных признаков объектов в одном классе не превышает eдоп(см. рис.7.2).

Значение величины eдоп задается в процентах и зависит от желательной степени подробности искомых характеристик. Затем выбирается значение параметра b, которое должно удовлетворять условию

Рис. 7.2. Аппроксимация функции S(e )

Значения By и Bkопределяются следующим образом. При вычислении значений Ti = j(ri) для элементов сети получили множества {Ti, i = 1,..., n; n - число УК} ; { Tj, j = 1,..., l; l- число КС}, из которых определяются размахи СРП для УК и КС: и соответственно

Рассмотрим пример. Пусть внутри подмножества путей длиной в 3 транзита П(3) нужно выделить классы. Расхождение результативных признаков внутри каждого класса должно быть не более 5% . Также пусть при значении b = 6% выделено [image] = 6 для УК и [image] = 7 для КС . Тогда имеем:

Следовательно, в качестве порога существенности связей нужно взять а = 0.965, при котором максимальное значение критерия (7.4) будет соответствовать оптимальному разбиению множества объектов П(3) на классы.

7.3.3. Алгоритм двухступенчатой декомпозиции сети

После того как определены способы назначения рангови выбора порога существенности связей , укажем последовательность действий для разбиения множества объектов на классы с помощью двухступенчатой декомпозиции .

1. На первой ступени проводится ранжировка элементов сети (узлови каналов) по способу, изложенному в разделе “Назначение рангов УК и КС сети”.

1.1. По заданному значению e с помощью (7.7) определяется значение разницы составляющих результативного признака b.

1.2. При заданных значениях нагрузки элементов сети и полученном значении b проверяется условие (7.5); при выполнении условия элементам назначается одинаковый ранг. В результате получаем значения рангов для УК и КС, которые определяют диапазон изменения компонентов вектора Z = (Z(1),...,Z(2d+1))

1.3. По признаку, характеризующему длину пути, множество П разбивается на подмножества П(d), d = 1,...,D , где D - максимальное число транзитных участков для путей piМ П .

2. Этапы второй ступени декомпозиции выполняются для каждого подмножества П(d), d = [image].

2.1. По заданному допустимому значению eдоп и по (7.6) находится значение порога существенных связей а.

2.2. С помощью (7.2) строится матрица близости S = п п Sij п п , i, j=1,...,Nd; Nd- размерность подмножества П(d).

2.3. Выполняется непосредственно сама классификация. Для максимизации критерия (3.2) используются две последовательные процедуры:

А. ОБЪЕДИНЕНИЕ. Эта процедура применяется к разбиениям R(k)={ [image]} (k = 0,1,...). Начальное разбиение k = 0 состоит из Nd одноэлементных классов.

От разбиения R(k)... переходим к разбиению R(k+1)с помощью объединения двух классов R(k) так, чтобы дать максимальное приращение показателю f(а,R). Для этого анализируются связи f( [image] ) между классами и объединяются такие [image] и [image], что суммарная связь максимальна и положительна. Процесс прекращается, как только среди сумм f( [image] ) не будет положительных.

Б. ПЕРЕМЕЩЕНИЕ. Этот процесс применяется, начиная с разбиения [image], полученного в результате предыдущей процедуры. Он состоит в последовательном перемещении каждого объекта в порядке нумерации в тот класс, при котором показатель (7.4) получает наибольшее положительное приращение. Если перемещение любого объекта в другие классы не увеличивает f(а,R), процесс заканчивается. Сходимость алгоритма гарантируется монотонным ростом показателя (7.4) .

7.3.4 . Имитационная модель виртуального канала

Имитационная модель виртуального канала позволяет проводить детальное исследование процесса передачи и обработки пакетов на любом выделенном пути (маршруте сети), имитировать любые режимы функционирования ВК.

7.3.4.1. Концептуальная модель виртуального канала

Основной задачей анализа модели ВК является определение ВВХ доставки пакетов по исследуемому пути при заданных характеристиках элементов пути, параметрах входящих потоков, алгоритмах обработки и дисциплинах обслуживания очередей.

ВК представляет собой коммутационный канал, обеспечивающий транспортировку пакетов между двумя портами сети, т. е. является некоторым маршрутом в сети (рис 7.3), состоящим из последовательности n УК и (n-1) КС, по которому осуществляется передача информации из узла источника УК1 в узел-адресат УК n.

Особенностью имитационной модели ВК является отображение потоков, циркулирующих по сети и влияющих на процесс прохождения пакетов по выделенному (моделируемому) пути.

Поток (УК1, УК n), характеристики процесса доставки пакетов которого нас интересуют, будем называть основным (полезным - ПП). Потоки, маршруты которых проходят через один или несколько транзитных участков выделенного пути, будем называть транзитными. Транзитные потоки (ТП ) вносят задержки в процесс доставки пакетов ПП, занимая ресурсы выделенного пути.

Рис. 7.3. Структура модели виртуального канала

Пакеты ТП поступают на выделенный путь со смежных узлов, не принадлежащих рассматриваемому пути. Пакеты ТП, поступающие на i-й узел, I = [image] , в соответствии с маршрутом покидают выделенный путь после обработки в j-м узле, j = [image] . При j = i пакеты сразу покидают выделенный путь после обработки в узле, на который они непосредственно поступили. При j < i пакеты транзитного потока следуют встречно пакетам ПП. При j > i пакеты ТП совпадают по направлению с пакетами ПП. Таким образом, если выделенный ВК содержит n узлов коммутации, то на каждый из них поступает не более n транзитных потоков, различающихся дальнейшим маршрутом следования.

Формат пакета в модели целесообразно сделать одинаковым для пакетов всех потоков и предусмотреть в формате поля для фиксирования: момента времени поступления пакета в выделенный путь (timestart) (момент генерации пакета); времени очередного "особого" момента (timimmed) относительно момента генерации; времени обслуживания пакета; параметра, по которому определяется время обслуживания; приоритета пакета; адресата пакетов (номера узла выделенного пути, после обработки в котором пакет покидает моделируемый ВК); состояния пакета: активный (true), когда пакет может обрабатываться или передаваться, иначе - пассивный (false).

В адресное поле пакетов ПП записывается идентификатор J = n+1, пакетам ТП, поступающим на i-й узел и покидающим выделенный путь после обработки в j-м узле, присваивается идентификатор J = j, j = [image] . После обработки пакета на каждом узле проверяется адресное поле пакета. Если j-й адрес соответствует номеру узла или меньше, то пакет покидает выделенный путь; если j больше номера узла, то пакет отдается передатчику для передачи на следующий узел. В конце пути после обработки пакетов в n-м УК, если идентификатор j > n, то пакет направляется в сумматор пакетов ПП, доставляемых адресату для последующей статистической обработки. Такой способ идентификации пакетов позволяет унифицировать процедуру их обработки на каждом узле независимо от типа потока, которому принадлежит данный пакет.

На выходе пакета ПП из пути вычисляется разность (timimmed - timestart), которая характеризует время транспортировки пакета по выделенному пути.

Каждый пакет после его обработки или передачи из активного (true) состояния переходит в пассивное (false), в котором сохраняется для дальнейшего сравнения состояния со следующими за ним пакетами. Для генерации ТП нужно знать среднюю интенсивность и распределения моментов поступления пакетов со смежных узлов, не принадлежащих выделенному пути.

Интенсивность потока, поступающего на любой узел сети, можно определить по формуле: , где l0-интенсивность потока, поступающего от абонентских пунктов, непосредственно подключенных к данному УК, ln- интенсивность потоков из смежных узлов; l - число смежных узлов; b = 1 - Pn коэффициент пропускания входного потока; Pn - вероятность потери сообщения на смежном узле из-за возможного переполнения буферного накопителя узла.

Если вероятность потерь в УК из-за переполнения его входного либо выходных буферов не превышает 10-9, то можно положить Pn =0 и определить интенсивность в виде:

. (7.8)

Моделируемый путь содержит n УК. На каждый узел поступает ТП с интенсивностью [image],определяемой по (7.8). В свою очередь каждый ТП является суперпозицией потоков с интенсивностями [image], где [image] - интенсивность потока, поступающего от абонентских пунктов, непосредственно подключенных к данному УК выделенного пути и покидающих путь после обработки в j-м УК; [image] - интенсивность потоков из смежных узлов, сообщения которых покидают выделенный путь после обработки в j-м узле; l-число смежных узлов с i-м УК выделенного пути; [image] - суммарная интенсивность ТП, поступающего в i-й УК выделенного пути, пакеты которого покидают путь после обработки в j-м УК.

Имея матрицу трафика сети и таблицы маршрутизации потоков, можно по (7.8) однозначно определить интенсивность входящего потока на любом узле сети и сформировать интенсивности транзитных потоков в соответствии с алгоритмом, приведенным далее.

Распределение моментов поступления пакетов ТП неизвестно, обычно на основе предельных теорем его принимают пуассоновским. Учитывая, что процесс транспортировки сообщений по сети управляется линейными и сквозным протоколами, можно полагать, что межузловым потокам присуща определенная регулярность. Необходимую степень регулярности можно отразить, задавая моменты поступления пакетов распределением Эрланга. Порядок Эрланга К для генерируемых ТП можно варьировать в широком диапазоне. Значение К = [image] соответствует регулярным ТП со случайно сдвинутыми моментами поступления пакетов на разных УК.

Буферные накопители в модели ВК размещаются как перед УК, так и перед КС. В случае, если при поступлении пакетов в УК или КС последний оказывается занятым, поступившие пакеты ставятся в очередь во входном буфере. Пакеты ТП после обработки на узле ВК, если их адрес соответствует данному узлу, покидают выделенный путь и больше не оказывают влияния на остальные пакеты.

Ко входу каждого УК ВК подключается генератор ТП, генерирующий поток с интенсивностью [image] ; [image] получаем в результате выполнения алгоритма формирования интенсивности ТП, который приводится далее. К первому узлу коммутации выделенного пути от отдельного генератора поступают также пакеты ПП.

Модель ВК строится по принципу особых состояний, которые необходимы для правильного воспроизведения времени доставки пакетов. В силу того, что УК и КС представлены одинаковой формализованной схемой, выделяется три особых состояния: возникновение требования на передачу пакетов (приход пакетов в систему); возникновение требования к ожиданию в очереди (УК или КС); поступление требования на обслуживание прибором (УК или КС). Первое состояние для КС возникает при окончании обслуживания пакетов в УК и является окончанием обработки на УК.

Поиск состояний производится по их последовательности в модели. Модули состояний (МС) получают управление от программы-диспетчера, в которой анализируется изменение каждого состояния и отражается в базе данных имитационной модели. После окончания работы модуля обработки в соответствии с последовательностью событий программа-диспетчер передает управление следующему модулю обработки состояния.

Таким образом, исходя из основной задачи анализа, связанной с оценкой ВВХ процесса доставки пакетов между корреспондирующими узлами, в имитационной модели ВК учитываются состояния очередей узлов и каналов, дисциплины обслуживания очередей, действия транзитных потоков, отображающих влияние на выделенный путь нерассматриваемой в модели части сети. Способность модели ВК настраиваться на любые направления обмена дает возможность использовать ее при построении алгоритмов маршрутизации, при параметризации конкретного пути обмена, при параметризации протоколов.

7.3.4.2. Алгоритм формирования интенсивностей транзитных потоков

Интенсивности для каждого из n ТП формируются поочередно согласно (7.8). При этом исходными данными являются:

1. Матрица тяготений L = ||l ij|| , i, j = [image] (N - число узлов), в которой заданы интенсивности потоков, поступающих в сеть. Нумерация потоков упорядочена согласно возрастанию номеров узлов-источников. Общее количество потоков при наличии обмена информацией каждого УК с каждым равняется N(N - 1).

2. Таблица маршрутизации потоков задана в виде матрицы размерности N(N - 1)(D + 1), в которой число строк равно числу потоков, число столбцов - числу узлов в максимальном по длине маршруте. В каждой строке матрицы записаны номера узлов, составляющих основной маршрут потока, номер которого соответствует номеру строки. Для потоков, маршрут которых состоит из числа узлов менее (D + 1), свободные столбцы заполнены нулями. Узлы, записанные в первый столбец матрицы маршрутизации, будем называть узлами, входящими в маршрут на первом уровне, узлы 2-го столбца - соответственно узлами 2-го уровня и т.д.

Пусть задан конкретный путь обмена информацией, состоящий из последовательности n узлов (n - 1 транзитов). Соответственно для i-го УК, i = 1,..., n нужно определить интенсивности (n- i+1) потоков: потока, проходящего через узел ni и поступающего встречно полезному потоку (эти потоки вносят задержки, поскольку занимают ресурсы выделенного пути на i-м узле перед их уходом с трассы выделенного ВК), потока, проходящего через узлы ni, (ni+1) и т.д., потока, проходящего через узлы (ni, (ni+1),..., nn).

Поток, поступающий на узел ni, покидает маршрут после nj узла с вероятностью

. (7.9)

Значения интенсивностей ТП, проходящих через узел ni рассматриваемого пути, и вероятности их ухода вычисляются одновременно по следующей схеме:

1. Формируется вектор суммарных интенсивностей потоков, поступающих на узел ni и покидающих узел nj, [image] ; идентификатором j = 0 обозначается поток, получаемый по (7.8), идентификатором j > 0 обозначаются суммарные потоки, поступающие на ni и покидающие nj.

2. Просматривается строка матрицы маршрутизации потоков.

3. Если узел ni входит в маршрут потока, то необходимо проверить является ли узел продолжением рассматриваемой строки матрицы маршрутизации потоков (ММП). Если да, то снова возвращаемся к п.2, в противном случае переходим к п.4 .

4. Если узел ni не входит в маршрут потока, соответствующего рассматриваемой строке, переходим к п.7.

5. В зависимости от того, какое звено заданного пути входит в рассматриваемый маршрут (звено ni, либо звено (ni, (ni+1)), либо звено (ni, (ni+1),..., nn), формируется интенсивность ТП, проходящего через это звено.

6. Если узел ni входит в рассматриваемый маршрут на j-м уровне, то интенсивность потока, соответствующего рассматриваемой строке ММП, суммируется с интенсивностью ТП, на которую указывает идентификатор "0", а также с интенсивностью, на которую указывает идентификатор j.

7. Переходим к следующей строке ММП и повторяем действия, начиная с п.2

По приведенному алгоритму определяются интенсивности ТП для каждого узла ni, [image] заданного пути обмена информацией и их вероятности ухода с пути, определяемые по (7.9).

7.3.4.3. Особенности программной реализации
модели виртуального канала

Моделирующая программа, как и всякая обычная процедура, является последовательно выполняемым набором операторов. Узлы в моделирующей программе ВК представляются тремя основными устройствами (рис.7.4): первое обеспечивает перенос пакетов между узлами (модуль Removpaket); второе обеспечивает коммутацию пакетов - это СМО с одним обслуживающим прибором, одной входной очередью и более чем одной выходной очередью (модуль УК); третье обеспечивает передачу пакетов следующему узлу - это СМО с одной входной очередью (модуль КС).

Наряду с разработкой алгоритмов EVT, обеспечивающих ускорение моделирования ВК методом Монте-Карло, ускорить процесс моделирования можно и за счет рационального построения программной реализации имитационной модели (ПРИМ) ВК.

Рис.7.4. Структура моделирующей программы ВК

Обычно в ПРИМ используется механизм системного времени, в основе которого лежит принцип особых моментов, - на очередном шаге выполняется то особое состояние, время наступления которого ближайшее к текущему моменту.

Для ускорения процесса моделирования в ПРИМ ВК понятие механизма системного времени можно свести к механизму системного перехода, который обеспечивает переход к нужному модулю и позволяет выполнять несколько особых состояний при очередном просмотре модулей ПРИМ. Например, момент поступления на передачу пакета в i-м узле предшествует моменту Ti+1 поступления в (i+1)-й узел транзитного пакета. Оба эти события можно выполнить одновременно, если момент завершения передачи в i-м узле более поздний, чем момент Ti+1 .

Ускорить ПРИМ возможно также за счет сокращения количества процессорных циклов при продвижении очередей. Например, за переносом пакета из буфера FIFO на обслуживание в УК (КС) следует продолжительный процесс продвижения всех пакетов, находящихся в буфере. Если средние длины очередей в буферах [image], то для обслуживания n узлов требуется (n· [image] ) дополнительных операций. Представление буфера в виде кольца (ящики пакетов связываются друг с другом через их адреса) позволяет при переносе пакета на обслуживание заменять адрес первого пакета в кольце, т.е. вместо подвижки всех пакетов в буфере адрес первого пакета заменяется на адрес второго пакета. Такой способ сокращает количество процессорных циклов, затрачиваемых на обработку одного пакета, в [image] раз на каждом узле.

Представленная структура позволяет выполнять нужные модули и обрабатывать несколько особых состояний при очередном просмотре модулей ПРИМ.

Модель состоит из программных модулей, реализующих алгоритмы обслуживания моделирования, и модулей, реализующих алгоритмы обработки имитируемых физических действий.

Алгоритмы обрабатывающих модулей поддерживают имитацию процесса движения пакетов по маршруту. Если обратиться к распространенной формализации поведения систем с параллельными и асинхронными процессами в виде сетей Петри, то функционирование обрабатывающего модуля может быть соотнесено со срабатыванием некоторого составного перехода.

Алгоритм Removpaket обеспечивает переход пакетов между узлами. В нем определяются особые состояния - возникновение требований к передаче пакетов. Функциональная схема алгоритма показана на рис. 7.5 Алгоритм выполняет следующие действия: Т i1 - вызывается генератор ТП (Genertranzit), и, если в интерфейсе генератора отсутствует пакет, генерируется транзитный пакет; Т i2 - просматривается генератор полезного потока, в случае отсутствия пакета в интерфейсе генератора генерируется пакет ПП; Т i3 , Тi5 - регистрируется транзитный пакет в интерфейсе Removpaket; Т i4 - регистрируется полезный пакет в интерфейсе Removpaket; Т i6 - в Removpaket регистрируется пакет с выхода КС предыдущего узла; Тд - управление передается в программу-диспетчер

Рис.7.5. Схема функционирования устройства Removpaket

Рис. 7.6. Функциональная схема устройства УК и КС

Условия имеют следующий смысл: Р i1=1 - если i=1, иначе i і 2; Рi2=1 - момент появления транзитного пакета предшествует появлению полезного пакета, иначе наоборот; Р i3=1 - есть пакет в КС предыдущего узла, иначе управление передается в программу-диспетчер; Р i4 - момент появления транзитного пакета предшествует моменту появления пакета из КС предыдущего узла, иначе наоборот.

Схема функционирования УК и КС показана на рис. 7.6. В ней использованы три основные процедуры: Waitpaket обеспечивает регистрацию пакетов в очереди УК (КС); Servicpacet обеспечивает обслуживание пакетов из очереди или непосредственно без ожидания; Servpriority обеспечивает обслуживание приоритетных пакетов.

Соответственно схеме: Р i1= 1 - момент поступления пакета в УК (КС) предшествует моменту окончания обслуживания предыдущего пакета в УК (КС), иначе наоборот; Р i2=1 - обслуженный пакет передан следующему устройству, иначе управление передается в программу-диспетчер; Р i3= 1 - буфер не пустой, иначе обслуживается поступающий пакет без ожидания; Р i4= 1 - модель работает с приоритетом, иначе модель работает без приоритетов; Т i1 - регистрируется поступление пакетов в очередь УК (КС); Т i2 - обслуживается поступающий пакет без ожидания; Т i3 - выбирается приоритетный пакет из буфера; Т i4 - обслуживается пакет; Т i5 - поднимается уровень приоритетов пакетов, стоящих в очереди перед выбранным в (Т i3) к обслуживанию; Тi6 - переход к обслуживанию первого пакета, стоящего в очереди; Тд - управление передается в программу-диспетчер.

Рассмотренные схемы функционирования межузлового устройства, УК, КС повторяются одинаково для всех узлов ВК с изменением идентификатора узла i.

К алгоритмам обслуживания моделирования относятся алгоритмы: построения генераторов транзитных потоков, механизма системного перехода, ведения базы данных, сбора и обработки статистик.

Механизм ведения системного перехода заключается в том, что программа-диспетчер просматривает активизации каждого устройства начиная с первого в первом узле и заканчивает последним в последнем узле. Под активизацией понимается занятие устройства пакетом. В таком случае идентификатор активизации ставится в положение (true). После завершения перехода пакета в следующее устройство идентификатор ставится в положение (false).

Для каждого узла предусмотрено три идентификатора активизации: первый - для межузлового устройства (Removpaket), указывает на появление пакета либо из предыдущего узла (для первого узла - из генератора полезного потока), либо из транзитного потока; второй - для устройства УК (Servicpaket-УК), указывает на занятость УК обслуживанием пакета; третий - для устройства КС (Servicpaket-КС), указывает на занятость КС передачей пакета.

Программа-диспетчер передает управление в очередное устройство при условии активизации предыдущего, в нем проверяет системное время, которое соответствует особым состояниям, присущим устройству, и передает управление особому состоянию, имеющему минимальное время поступления.

Removpaket следит за генерацией транзитных пакетов, сравнивая моменты их поступления с моментом поступления пакета из предыдущего узла (для первого узла - поступление пакета из генератора ПП) и регистрирует пакет, имеющий минимальное время поступления, после чего восстанавливает его активизацию, сообщая о требовании к передаче пакета.

Servicpaket и Waitpaket следят за регистрацией поступающих пакетов в очереди устройства УК и КС или за их обслуживанием в УК и КС, сравнивая момент поступления очередного пакета и момент окончания обслуживания предыдущего.

Программа-диспетчер следит за удалением транзитных пакетов, непрекращающих продвижение по маршруту после их обслуживания в УК.

Ведение баз данных включает следующие основные операции над данными:

- регистрацию генерируемого или поступающего в межузловое устройство пакета (Removpaket), возобновление активизации устройства;
- перенос пакета из межузлового устройства в очередь устройства УК (ОУК) и деактивизацию межузлового устройства;
- выбор пакета из ОУК для обслуживания и активизации устройства УК; в модели с приоритетом повышаются приоритеты пакетов, стоящих перед выбранным;
- перенос пакета с устройства УК в очередь КС (ОКС), если адрес пакета направляет его к следующему узлу; деактивизацию устройства УК;
- выбор пакета из ОКС для передачи и активизации устройства КС;
- перенос пакета из КС в устройство межузлового перехода и деактивизацию устройства КС.

Сбор и обработка выборочных данных выполняются процедурами, которые дополняют друг друга для получения основных характеристик ВК. Процедуры обеспечивают:

- сбор информации для вычисления вероятностей потерь в каждом буфере, входящем в ВК (в УК и КС), построение распределения потерь в зависимости от емкости буфера. В этих процедурах использована методика, предложенная в р. 7.3;
- получение необходимой информации при применении метода EVT для построении распределения вероятностей потерь в буфере конечной емкости по (3.15), (3.16);
- сбор времен доставки пакетов и их хранение для построения характеристик времени транспортировки пакетов по ВК. Выходными продуктами этих процедур являются: среднее время пребывания (транспортировки) пакетов в ВК; значение времени, за которое с вероятностью 0.99 пакет будет передан по ВК (99% времен пребывания пакетов в ВК не превышают определенного значения), и распределение времени пребывания (транспортировки) пакетов в ВК.

К оглавлению

7.4. Модель и методика расчета малых вероятностей потерь
в буфере конечной емкости

Аналитические методы расчета вероятности потерь, обусловленных переполнением буфера конечной емкости (далее - конечный буфер ( КБ)), основаны на модели очереди в СМО с буфером неограниченной емкости (далее - бесконечный буфер (ББ)), и вероятность потерь рассчитывается, исходя из функции распределения очереди F(h) в ББ в виде Р = 1- F(h Ј N), где N - емкость КБ. Однако такая модель не адекватно отражает процесс потерь в КБ, что приводит с возрастанием нагрузки к значительным погрешностям при оценке вероятности потерь в виде Р = 1- F(h) [ГИК]. Излагаемый метод расчета малых и очень малых (Ј 10-9) значений вероятности потерь в КБ опирается на модель поведения очереди в КБ, которая позволяет построить статистику и получить выражение для оценки вероятности потерь. При этом время имитации сокращается в 4- 5 декад по сравнению с классическим методом Монте-Карло [Итал].

7.4.1. Модель процесса потерь в буферном накопителе
конечной емкости

Анализ процесса потери заявок при переполнении КБ СМО можно свести к анализу наращивания очереди в буфере и превышению ею емкости КБ, что и приводит к потерям.

В стационарной СМО (коэффициент загрузки r < 1) буфер время от времени освобождается, т.е. длина очереди время от времени становится равной нулю. В моменты освобождения можно считать, что система находится в первоначальном состоянии и промежуток занятости (период регенерации) есть промежуток между двумя моментами освобождения.

События, происходящие в системе в разные промежутки занятости, независимы и подчиняются одним и тем же закономерностям. Из этого следует, что в силу стационарности СМО анализ происходящих событий в одном промежутке занятости эквивалентен для всех промежутков.

В течение интервала занятости заявки поступают и покидают систему. Если во время занятости обслуживающего прибора в систему поступает более одной заявки (условие 1), то очередь в буфере начинает наращиваться. И, наоборот, очередь сокращается, если за время обслуживания в систему поступает не более одной заявки (условие 2). Сокращение очереди продолжается до тех пор, пока не наступает условие 1, либо пока не освобождается полностью буфер и кончается промежуток занятости.

Количество процессов наращивания и процессов сокращения очереди, происходящих в одном промежутке занятости, одинаково.

Для оценки величины потерь в КБ рассмотрим два буфера - один с конечной емкостью N и другой с бесконечной (неограниченной) емкостью - ББ.

Пусть имеются наблюдения событий, происходящих в ББ и КБ, при одной и той же траектории поступления и обслуживания заявок.

Введем следующие обозначения: N - емкость КБ; L - максимальное значение длины очереди в ББ после i-го процесса наращивания; [image] [image] - значения длины очереди после i-го сокращения в ББ и КБ соответственно (перед (i+1)-м наращиванием); z2, [image] - количество поступивших в систему заявок, обеспечивающих i-е наращивание; Wi, [image] - количество обслуженных заявок при i-м сокращении; Пi, [image] - количество потерянных заявок в КБ за счет i-го наращивания; П(n) - количество потерянных заявок при n наращиваниях на интервале занятости (регенерации) .

Развитие событий в периоде регенерации имеет два различных варианта:

первый - все максимальные значения длины очереди в ББ Li, [image], n і 1, не превышают емкости КБ, т.е. LiЈ N "i. В этом случае потерь заявок в КБ нет, события как в КБ, так и в ББ развиваются аналогичным образом;

второй - несколько максимальных значений длины очереди, либо все значения в ББ [image] , n і 1, оказались больше, чем емкость КБ. Для определения потерь заявок в этом случае рассмотрим ряд ситуаций (рис. 7.6 - 7.8).

На рисунках можно выделить два варианта: 1-й - когда [image] 2-й - когда [image], что соответствует ситуации, когда устройство обработки ждет поступления заявок для обслуживания.

Рассмотрим ряд ситуаций при

1. n=1. Если в промежутке произошел только один процесс наращивания до значения L1(рис. 7.6), то потери в КБ составят (L1- N) заявок, т.е.

Рис. 7.7. Однократное превышение емкости N на интервале регенерации

2. n=2. На рассматриваемом периоде занятости кроме первого процесса наращивания L1>N возникает еще один процесс наращивания очереди к L2. Случай L2Ј N к дополнительным потерям в КБ не приводит и с точки зрения расчета потерь в конечном буфере ничем не отличается от вышерассмотренного (n = 1) .

Интерес представляет случай L2>N. Рассмотрим реализацию L2>N, L1> L2 (рис. 7.8).

Рис. 7.8. Двукратное превышением емкости N
на интервале регенерации L1> L2

После первого наращивания в ББ процесс сокращения освободит буфер до уровня [image], а с момента начала процесса второго наращивания до уровня L2 в систему поступит z2 заявок. Такое же количество заявок поступит в КБ. Однако в КБ процесс второго наращивания начинается с уровня [image]. Количество потерь П1 за счет первого наращивания к L1 есть П1= ( L1- N) . Второе наращивание до уровня L2 происходит за счет поступления в систему [image] заявок. В КБ поступление этого же количества заявок за счет второго процесса наращивания приводит к накоплению заявок до уровня, меньшего N: , поскольку L2 < L1.

Таким образом, при рассматриваемой реализации (L1, L2 > N, L1 > L2) процесса второго наращивания дополнительных потерь в КБ не происходит, т.е. П2 = 0;L2 > L1.

Пусть теперь при втором наращивании имеет место L1, L2 > N, L2 > L1. Рассуждаем так же, как при L2 > L1. За счет второго процесса наращивания в КБ будет достигнут уровень [image] поскольку L2 > L1. При такой реализации за счет процесса второго наращивания в КБ возникают потери в количестве П2= N + (L2 - L1) - N = L2 - L1 заявок, т.е.

Таким образом, общее количество потерянных заявок составит [image] или П1+ П2 = L1 - N + L2 - L1 - L2 - N; L2 > L1.

Следовательно, и при этой реализации развития процессов наращивания число потерянных заявок в КБ совпадает с максимальным превышением уровня N в бесконечном буфере, т.е.

3. n=3. Пусть при третьем процессе наращивания в ББ достигнут уровень L3, превышающий значение N, т.е. L1,L2,L3 > N (рис.7.9).

Рассуждая аналогично случаю двухкратного наращивания, можно показать, что в КБ количество потерянных заявок П(n=3) = П(n=2) + П3, где П3- количество потерянных заявок за счет третьего наращивания. Из предыдущего выделим два случая - L 1і L 2 и L 1< L 2:

а) L1і L2 Ю П(n=2) = L1 - N можно показать, что в КБ третий процесс наращивания до уровня L3 начинается с уровня [image] , где [image] и количество потерянных заявок за счет третьего наращивания П3 будет составлять [image], если L3 > L1 и П3= 0 , если L3 < L1.

Рис. 7.9. Трехкратное превышение емкости N
на интервале регенерации, L3>L1>L2

Общее количество потерянных заявок

б) L1 < L2 Ю П(n=2) = L2 - N в этом случае третий процесс наращивания до уровня L3начинается с уровня [image] и количество потерянных заявок за счет третьего наращивания П3составит [image], если L3 < L2.

Общее количество потерянных заявок

Следовательно, количество потерянных заявок в КБ при трехкратном превышении емкости N будет составлять

Аналогичные рассуждения показывают, что и при n- кратном процессе наращивания на одном промежутке занятости потери в КБ будут определяться максимальным превышением значения N емкости в ББ на одном промежутке занятости системы.

Случай [image], когда траектории развития процессов в КБ и ББ расходятся, возникает исключительно редко, и его влиянием при расчете вероятности потерь PП можно без ущерба для точности конечного результата пренебречь.

Таким образом, из анализа модели СМО с ББ следует, что потери в СМО с КБ определяются максимальным значением очереди, достигаемым в ББ на одном промежутке занятости независимо от траекторий развития процессов наращивания и сокращения.

И, следовательно, количество теряющихся в КБ заявок из-за переполнения его емкости равно разности между максимальной величиной очереди, образуемой в ББ на промежутке занятости, и значением емкости КБ.

Максимальная величина очереди в ББ, достигаемая в процессе наращивания, не зависит ни от количества процессов наращивания и сокращения, ни от количества заявок, обслуженных в этом промежутке, т.е. не зависит от длины промежутка занятости. Обозначим через Li максимальное значение очереди в ББ, достигнутое на i-м промежутке занятости. Значение Liне зависит от номера промежутка занятости, а зависит только от временных интервалов между поступлениями и интервалов обслуживания. Поскольку эти времена являются случайными величинами, то Li" i тоже случайная величина. При существовании бесконечного количества промежутков занятости можно утверждать, что величина L имеет функцию распределения вероятности своих значений.

Теперь можно сказать, что если в течение промежутка занятости происходит потеря одной заявки в КБ емкостью N, то это соответствует тому, что в ББ максимальная величина очереди в этом промежутке достигает значения L = N + 1; если теряется две заявки, то L = N + 2 и т.д.

Таким образом, можно утверждать, что вероятность потери заявок на конечном буфере емкостью N есть функция от хвоста распределения значений максимальной величины L, начиная со значения L=N+1. Следовательно, если функция распределения при L=N есть F(x), то вероятность потерь есть

Pn(N) = f[Q(N + 2)] = f[1 - F(N)].

7.4.2. Методика расчета вероятности потерь в КБ

Изложим основные шаги методики сбора и обработки выборочных данных при построении модели для оценивания потерь в КБ методом экстремальных статистик (EVT) [Galambosch?Muf]. Полагаем, что время моделирования или объем общей выборки и признак выделения точек регенерации определены на более ранних стадиях моделирования.

Итак, в процессе сбора и обработки данных необходимо:

  1. Выделить интервалы регенерации из общей выборки значений очереди, подсчитать число интервалов и оценить их среднюю длину [image] ( [image] - среднее число заявок, поступивших в систему на интервале регенерации. При поступлении каждой заявки фиксируется значение очереди в системе).
  2. По всем интервалам выделить максимальные выборочные значения, которые образуют последовательность независимых значений случайной величины - длины возможной очереди в КБ - для последующей обработки методом EVT.
  3. Образовать порядковую последовательность из независимых значений длин очереди и оценить вероятности этих значений.
  4. Разделить полученную в п.2 последовательность на M подпоследовательностей, содержащих по n выборочных значений (n можно принимать равным 50 значениям).
  5. Образовать последовательность из максимальных значений xi i=1,.., М подпоследовательностей.
  6. Образовать порядковую последовательность из максимальных значений (п.5), оценить вероятность каждого значения в ней и построить в виде гистограммы оценку [image] функции распределения максимального значения выборочной последовательности (п.2):
    где hi есть количество повторов i-го выборочного значения.
  7. Образовать последовательность значений для каждого xi.
  8. Оценить значение an как статистику по формуле
     ,
    где [image] , [image], N - число шагов гистограммы, и значение bn по выражению [image] .
  9. Исходя из допустимой вероятности потерь, определить необходимую емкость N КБ из выражения
    PN     = [image] [image],
    где 2Dx = hi - hi-1 - длительность интервала между значениями hi, hi-1 дискретной величины.
  10. Если допустимая вероятность задана не малым значением, то необходимую емкость буфера можно определить, исходя из выражения:

PN =

7.4.3. Замечание к определению независимых периодов регенерации

Существует много представляющих интерес случаев регенерирующих процессов, в которых среднее время между последовательными точками регенерации конечно, но, возможно, недостаточно мало, чтобы в ходе моделирования можно было пронаблюдать разумное число моментов регенерации.

В этих случаях можно использовать следующий подход к регенерации для моделирования процессов с большим временем возвращения в состояние S.

1. Выбрать такую (малую) область А, содержащую S, что среднее время возвращения в А конечно.

2. Предположить, что процесс приближенно регенерируется каждый раз при возвращении в А.

На вопрос о том, в какой степени может быть свободен выбор области квазирегенерации А, чтобы хорошее приближение все еще сохранялось, определенного ответа нет. Ограниченный опыт применения этого метода дает основание предположить, что приближения удовлетворительны при довольно-таки свободном выборе таких областей.

Экспериментально на имитационных моделях СМО двух типов MЅ DЅ 1 и MЅ MЅ 1 установлено, что при оценке вероятности потерь по изложенному способу можно определять точки регенерации не только как S={0}, но и S={0,1}, S={0,1,2} и т.д. Главное, чтобы сохранялась независимость событий, происходящих в следующем периоде от событий в предыдущем периоде. Тогда полностью сохраняется справедливость модели потерь в КБ, основанной на распределении максимальных значений длины очереди на периодах регенерации в ББ. При нагрузках r і 0.5 время имитации для S={0,1,2} существенно уменьшается по сравнению с S={0}при получении одного и того же количества выборочных значений, используемых в расчете параметров EVT.

К оглавлению


МОДЕЛИРОВАНИЕ ТЕЛЕКОММУНИКАЦИОННЫХ СЕТЕЙ

Министерство Российской Федерации по связи и информатизации
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ТЕЛЕКОММУНИКАЦИЙ им. проф. М.А. БОНЧ-БРУЕВИЧА

О.И. Кутузов, Т.М. Татарникова

МОДЕЛИРОВАНИЕ
ТЕЛЕКОММУНИКАЦИОННЫХ СЕТЕЙ

Учебное пособие

Изложен метод статистического моделирования, включающий сведения о датчиках случайных чисел, моделирования случайных факторов, планирования статистического эксперимента, и аналитико-статистический подход к имитационному моделированию сложных систем, какими являются телекоммуникационные сети. Приведены примеры моделей для оценивания вероятностно-временных характеристик сетей. Предназначено для студентов, аспирантов и инженеров.

СОДЕРЖАНИЕ :

Предисловие

1. Концепция статистического моделирования
1.1. Технологическая схема имитационного моделирования систем
1.2. Концепция статистического моделирования

2. Построение датчиков БСВ
2.1. Датчики БСВ
2.2. Метод середины квадрата
2.3. Мультипликативный конгруэнтный метод
2.4. Тестирование равномерности
2.5. Тестирование независимости

3. Моделирование случайных факторов
3.1. Случайные факторы и их модели
3.2. Моделирование случайных событий
3.3. Моделирование дискретных случайных величин
3.4. Моделирование непрерывных с.в. методом обращения
3.5. Моделирование экспоненциальной с.в.
3.6. Доказательство метода обращения
3.7. Линейные преобразования с.в.
3.8. Моделирование равномерной с.в.
3.9. Моделирование эрланговской с.в.
3.10. Моделирование нормальной с.в.

4. Планирование статического эксперимента
4.1. Задача планирования
 4.2. Теоретическое решение задачи планирования
 4.3. Метод автоостанова
4.4. Интерактивный контроль

5. Статический расчет надежности
5.1. Классическая задача надежности
 5.2. Аналитический метод расчета
5.3. Статистический метод расчета
5.4. Пример статистического расчета
5.5. Сравнение аналитического и статистического методов

6. Аналитико-статистический подход
6.1. Общая характеристика
6.2. Метод взвешивания
6.3. Пример взвешенного моделирования
6.4. Метод расслоения
6.5. Примеры расслоения эксперимента

7. Реализации аналитико-статистического подхода применительно к анализу телекоммуникационных сетей
7.1. Задачи системного проектирования ТКС
7.2. Статистическое оценивание функциональной надежности сети
7. 3. Аналитико-статистический расчет вероятностно-временных характеристик транспортировки пакетов по сети
7.4. Модель и методика расчета малых вероятностей потерь в буфере конечной емкости

Заключение

Литература


ЛИТЕРАТУРА

Назад

Содержание

ЛИТЕРАТУРА

  1. Башарин Г.П. Модели информационно-вычислительных систем. М.: Наука, 1993. – 156 с.
  2. Бухштаб А.А. Теория чисел. – М.: Просвещение, 1966. –384 с.
  3. Вероятностные методы в вычислительной технике: Учеб. пособие для вузов по спец. ЭВМ /Под ред. А.Н.Лебедева и Е.А.Чернявского. – М.: Высш. школа, 1986. – 312 с.
  4. Гаскаров Д.В., Истомин Е.П., Кутузов О.И. Сетевые модели распределенных систем – СПб.: Энергоатомиздат, Санкт-Петербургское отделение, 1998 – 353 с.
  5. Голованов О.В. и др. Моделирование сложных дискретных систем на ЭВМ третьего поколения. Опыт применения GPSS. – М.: Энергия, 1978. – 178 с.
  6. Ермаков С.М., Михайлов Г.А. Статистическое моделирование. - 2-е изд., доп. – М.: Наука, 1982. – 296 с.
  7. ЕС ЭВМ. ППП для моделирования и исследования на ЭВМ дискретных систем (GPSS/360). Руководство программиста ПРО 309.007 Д2. 1977. – 51 с.
  8. Журбенко И.Г., Кожевникова И.А., Клиндоухова О.В. Определение критической длины последовательности псевдослучайных чисел // Вероятностно-статистические методы исследования. – М.: Изд-во МГУ, 1983. – С. 18–39.
  9. Задорожный В.Н. "Статистическое моделирование": Уч. пособие. – Омск, 1996. – 92 с.
  10. Иглхарт Д.Л., Шедлер Д.С. Регенеративное моделирование сетей массового обслуживания / Пер. с англ.; Под ред. В.В.Калашникова. – М.: Радио и связь, 1984. – 136 с.
  11. Исследование операций / Пер. с англ.; Под ред. Дж. Моудера и С. Эл маграби. - М.: Мир, 1981. Т.1 – Методологические основы и математические методы. – 712 с. Т.2 – Модели и применения. – 677 с.
  12. Клейнен Дж. Статистические методы в имитационном моделировании/Пер с англ.; Под ред. Ю.П.Адлера и В.Н.Варыгина. - М.: Статистика, 1978. – Вып.1 – 221 с.; Bып.2 – 335 с.
  13. Клейнрок Л. Теория массового обслуживания / Пер. с англ.; Пер. И.И.Грушко; ред. В.И.Нейман. – М.: Машиностроение, 1979. – 432 с.
  14. Клейнрок Л. Вычислительные системы с очередями / Пер. с англ.; Под ред. Б.С.Цыбакова. – М.: Мир, 1979. – 600 с.
  15. Кнут Д. Искусство программирования для ЭВМ.Т 2. Получисленные алгоритмы. – М.: Мир, 1977. –728 с.
  16. Коваленко И.Н., Кузнецов Н.Ю. Методы расчета высоконадежных систем. – М.: Радио и связь, 1988. – 176 с.
  17. Кузнецов Н.Ю. Аналитико-статистический метод нахождения нестационарного коэффициента готовности у альтернирующего процесса восстановления //Кибернетика. – 1986. – N1. – С. 125–127.
  18. Кутузов О.И., Задорожный В.Н. Аналитико-статистический метод для расчета высоконадежных систем связи // Техника средств связи. Техника проводной связи. –1990. – Вып. 1. – С. 121–130.
  19. Нейлор Т. Машинные имитационные эксперименты с моделями экономических систем. – М.: Мир, 1975.
  20. Плакс Б.И. Расчет надежности систем со сложной структурой ускоренным методом Монте-Карло //Изв. АН СССР. Техническая кибернетика. – 1983. – N6. – С. 158–162.
  21. Полляк Ю.Г. Вероятностное моделирование на ЭВМ. – М.: Сов. радио, 1971. – 400 с.
  22. Потапов В.И., Задорожный В.Н., Кузнецова Е.М. Методология имитационного моделирования вычислительных систем: Учеб. пособие. – Омск, 1983. – 80 с.
  23. Потапов В.И., Задорожный В.Н., Кузнецова Е.М. Программные средства моделирования вычислительных систем: Учеб. пособие. – Омск, 1984. – 80 с.
  24. Разработка системы моделирования дискретных процессов "ППП СМДП-В" (GPSS-V). Руководство программиста. Ч.1 – 101 с., Ч.2 – 152 с., Ч.3 – 164 с., Ч.4 – 139 с., Ч.5 – 174 с.
  25. Рябинин И.Д., Черкесов Г.Н. Логико-вероятностные методы исследования надежности структурно-сложных систем.-М.:Радио и связь,1981. – 264 с.
  26. Соболь И.М. Метод Монте-Карло. – М., 1978. – 64 с.
  27. Соболь И.М. Численные методы Монте-Карло. – М., 1973. –212с.
  28. Шеннон Р. Имитационное моделирование систем – искусство и наука. – М.: Мир, 1978. – 418 с.
  29. Шрайбер Т.Дж. Моделирование на GPSS / Пер. с англ.; Пер. В.И.Гаргера, И.Л.Шмуйловича; Ред. М.А.Файнберг. – М.: Машиностроение, 1980. – 592 с.,
  30. Харин Ю.С., Степанова М.Д. Практикум на ЭВМ по математической статистике: Для мат. спец. ун-тов. – Минск, 1987. – 304 с.


ПРЕДИСЛОВИЕ

Назад

Содержание

Вперед

ПРЕДИСЛОВИЕ

Содержание пособия отражает разделы курса "Распределенные системы управления" (РСУ), в которых излагаются методы и модели статистического моделирования, используемого при анализе надежности и производительности телекоммуникационных систем. Ядром РСУ, условием выполнения в РСУ сетевой технологии являются телекоммуникационные сети (ТКС). Создание ТКС, организация и построение корпоративных и глобальных сетей требуют знания и активного применения средств анализа альтернативных решений.

Технология моделирования сложных технических систем, к классу которых относятся ТКС, опирается на имитационное машинное моделирование. Термин “имитационное моделирование” часто отождествляется с термином “статистическое моделирование”. Статистическое моделирование как метод машинной реализации имитационных моделей систем, подверженных случайным воздействиям, в настоящее время наиболее употребим.

Статистическое моделирование – это технология, основанная на применении законов математической статистики и способностях современных компьютеров порождать и обрабатывать за короткие промежутки времени огромное количество данных. Тем самым становится возможным строить достаточно точные модели реальных сложных объектов.

В пособии излагается материал, связанный со статистическим моделированием телекоммуникационных сетей, формально представляемых в виде случайных графов. Такой формализм широко используется в качестве математических моделей ТКС. Основное внимание уделяется двум аспектам статистического моделирования ТКС: способам генерации случайных факторов, определяющих состояния системы, и построению аналитико-статистических моделей для оценивания структурной, функциональной и временных характеристик сетей.

Содержание пособия несколько шире его названия. Фактически излагается материал статистического моделирования стохастических сетевых моделей. Последние являются распространенными формальными описаниями мультипроцессорных систем, транспортных систем, распределенных автоматизированных систем.


Заключение

Назад

Содержание

Вперед

ЗАКЛЮЧЕНИЕ

Статистическое моделирование позволяет находить решения там, где другие математические методы отступают в очевидном бессилии. С помощью статистического моделирования можно получать решения близкие к оптимальным для систем, с которыми реально приходиться иметь дело инженерам, управленцам, экономистам и многим другим категориям профессионалов.

Применение аналитико-статистических методов в задачах проектирования и эксплуатации современных сетей и подобных сложных систем обусловлена необходимостью определения дифференциальных и интегральных критериев качества функционирования сетей, в которых ряд ограничений задан алгоритмически с помощью имитационной модели, что требует больших затрат машинных ресурсов. Рассмотренный в пособии комплекс методов, моделей и алгоритмов аналитико-статистического моделирования в определенной степени способствует решению этой проблемы.

Предложенный подход к декомпозиции сети сохраняет взаимообусловленность параметров и процессов всей сети при имитационном моделировании конкретного пути обмена информацией. Прикладное значение данного результата состоит в методике ускоренного машинного анализа сети. Модель процесса потерь в буфере конечной емкости позволяет построить распределение потерь в КБ при произвольных законах поступления и обслуживания требований в системе. Имитационная модель доставки вызова по сети с альтернативными маршрутами отличается возможностью передачи повторных вызовов как от узла-источника, так и от промежуточных предыдущих узлов с альтернативными исходящими путями. Модель пригодна для анализа любого встречающегося в практических задачах графа сети. Имитационная модель ВК обеспечивает универсальность задания режимов работы ВК и реализацию метода максимальных статистик для расчета временных характеристик доставки пакетов. Аналитико-статистические методы обеспечивают расчет оценок вероятности потерь в буфере конечной емкости; распределения времени доставки пакетов по виртуальному каналу; вероятности связности (структурной надежности) сети; вероятности установления соединения на сети за время, не превышающее допустимое (оценка функциональной надежности сети); дифференциальных и интегральных временных характеристик доставки пакетов по сети.

Необходимо, однако, обратить внимание читателя на то, что любые методы исследования, в том числе и методы статистического и аналитико-статистического моделирования, требуют от нас большой работы и многократного активного контакта с предметом исследования. Подкрепим эту мысль словами самого Томаса Шрайбера,описавшего датчик БСВ, встроенной в язык GPSS: “Дополнительно следует указать, что после того, как модель создана, работа исследователя только начинается” [9].